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FOREWORD 


The  author  wishes  to  acknowledge  the  assistance  of  Mr.  Carl  S. 
Christensen,  who  did  much  of  the  original  logic  and  programing  of  SPURT. 
Mr.  Christensen,  now  with  the  Aerospace  Corporation,  Los  Angeles, 
California,  was  formerly  a  lieutenant  at  the  Air  Force  Special  Weapons 
Center. 
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ABSTRAC  T 


SPURT  is  a  five-degree>of-freedom  trajectory  digital  computer  program 
for  spinning  ui^fuided  space  probe  vehicles.  The  program  was  written  for 
the  Control  Data  Corporation,  1604  digital  computer. 

SPURT  will  compute  trajectories  for  a  vehicle  up  to  a  maximum  of  ten 
stages  and  has  provision  for  computing  the  trajectories  of  the  separated 
stages. 

This  generalised  program  computes  the  trajectory  over  an  oblate 
spheroidal,  rotating  Earth  with  atmosphere,  in  a  geocentric  rectangular 
coordinate  system.  All  input  and  output  data  are  in  geodetic  coordinates. 

Coasting  flight  trajectories  are  computed  in  two  subroutines.  The  first 
is  a  Keplerian  solution,  which  also  computes  orbital  elements  and  "look 
angles"  for  various  tracking  stations.  The  second  uses  three- degree- of- 
freedom  point  mass  equations  solved  by  numerical  integration. 

The  program  will  prepare  two  special  output  tapes.  One  is  used  in 
plotting  output  data  and  the  other  is  used  to  prepare  a  special  tape  for  the 
Atlantic  Missile  Range. 
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1.  INTRODUCTION. 

The  Spinning  Ungulded  Rocket  Trajectory  (SPURT'  digital  coniputer 
program  is  designed  to  provide  a  generalized  computer  program  for  calcula¬ 
ting  trajectories  of  spinning  unguided  space  probe  vehicles  such  as  the 
SLV-IB.  Many  trajectory  programs  are  available^*®  but  none  meet  the 
demands  required  by  the  Space  Vehicle  branch  (SWTTS)  of  AFSWC,  The 
program  is  written  in  a  mixture  of  FORTRAN  and  CODAP  for  use  on  the 
CDC-1604  computer  and  utilizes  subroutines  of  the  CO-OP  library.  The 
program  is  designed  to  handle  up  to  a  ten  '^*.ige  vehicle  for  both  powered 
and  unpowered  flight  with  provisions  to  calculate  the  trajectories  of  the 
separated  "expended"  stages. 

The  main  program  computes  the  powered  portion  of  the  trajectory  and 
controls  entry  into  the  various  subroutines.  All  unpowered  flight  trajectories 
are  computed  in  the  two  subroutines  "TWO  BOD"  and  "IMPACT."  Both  of 
these  subroutines  have  the  capability  to  calculate  the  trajectory  of  all  the 
separated  stages  and  of  the  payload.  TWO  BOD  is  a  Keplerian  trajectory 
program  with  provisions  for  calculating  look  angles  of  various  tracking 
stations,  while  IMPACT  integrates  the  equation  of  a  point  mass  with  drag. 


The  main  program  uses  a  flve-degree-of-freedom,  three-dimensional 
system  with  the  sixth  degree  constrained  by  a  table  of  spin  rates  that  are 
read  into  the  program.  Ani  oblate  rotating  earth  and  associated  gravitational 
potential,  standard  1962  atmosphere,  and  altitude- dependent  wind  provisions 
are  also  incorporated  into  the  program.  The  position  vector  is  calculated  in 
a  rotating  Earth-centered  coordinate  system,  while  the  angular  positions 
arc  calculated  in  a  launch-centered  coordinate  system. 
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The  procedure  includes  provisions  for  coasting  periods  between  stages, 
which  are  terminated  by  time.  Thrust  is  computed  from  thrust  vs.  time 
tables  and  corrected  for  atmosphere  back  pressure. 

Aerodynamic  forces  and  moments  about  the  center  of  mass  are  inter¬ 
polated  from  a  table  of  Mach  number  dependent  coefficients. 

Computations  are  carried  out  using  either  the  Adams  or  the  Runge- 
Kutta  method  of  numerical  integration.  Both  methods  can  be  used  in  either 
fixed  or  variable  step  size-mode. 

The  program  has  the  provision  for  writing  two  special  output  tapes. 

One,  a  plot  tape,  is  designed  to  be  used  with  a  special  plot  program  for  use 
on  the  AFSWC  plotter  and  can  plot  any  output  variable  against  any  other  output 
variable.  The  second,  a  BATT  tape,  is  used  with  the  BATT  program  to  pre¬ 
pare  magnetic  tapes  to  meet  specified  Atlantic  Missile  Range  formats. 

The  program  running  time  is  approximately  3  to  4  minutes  for  a  four- 
stage  vehicle  similar  to  the  SLV-lB.  The  integration  and  atmosphere 
routines  are  compiled  separately  in  octal  locations  6000  to  7514.  The  rest 
of  the  program  is  compiled  in  fixed  binary  mode  starting  at  octal  location 
10050.  The  two  parts  are  then  put  together  to  be  read  into  the  computer. 

This  method  has  the  advantage  that  the  integration  and  atmosphere  routines 
are  not  recompiled  every  time  a  change  is  made  in  the  rest  of  the  program. 

The  program  is  stored  in  core  according  to  the  following  octal  addresses. 


6000 

- 

7142 

Integration  Routine 

7300 

- 

7514 

Atmosphere  Routine 

7660 

- 

10026 

SQRTF,  EXPF,  &  LOGF  Routines 

10050 

- 

34102 

SPURT  Routine  (main  program) 

34103 

- 

34262 

SETTAB  Routine 

34263 

- 

34350 

EC  LOCK  Routine 

34351 

- 

34364 

SC  LOCK  Routine 

34365 

- 

56750 

TWO  BOD  Routine 

f-  subroutines 

56751 

- 

61516 

IMPACT  Routine 

61517 

- 

61615 

GEODED  Routine 

61616 

- 

62006 

ROTATE  Routine 

62013 

- 

67420 

LIBRARY  Routines 

70037 

_ 

71343 

COMMON 
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2.  EQUATIONS, 


SYMBOLS 


Dlmen«lon« 


A 

Axial  xnomont  of  Inartlal 

ft^  -  slug 

EXIT  araa  of  noaala 

In^ 

Aalmuth  aagla*  of  the  wind 

deg 

Ax 

Aalmuth  angle  of  the  axle 

deg 

•e 

EqMtortal  radiue  of  -the  Earth 

ft 

B 

Longitudinal  moment  of  inertia 

ft^-slug 

c4 

l/^yT-(2f  -  f^lain^  Sq  <Ref  6) 

ft 

Dreg  coefficient 

None 

^DB 

Powered  flight  drag  coefficient 

None 

^DC 

Coasting  flight  drag  coefficient 

None 

.N^  Normal  force  coefficient  with 

V  da  /  respect  to  angle  of  attack 

per  radian 

Moment  coefficient  with  respect  to 
angle  of  attack 

None 

Center  of  pressure  of  tho  missile** 

ft 

D 

Diameter  of  the  missile 

ft 

D 

Total  drag  on  the  missile 

d 

Derivative  of  a  variable 

None 

d 

Reference  length 

ft 

•e 

Eccentricity  of  the  Earth 

None 

? 

Total  force  vector  on  the  missile 

f 

Flattening  of  the  Earth 

None 

*  All  aslmttth  uifl«s  ar*  maasuvtd  clockwia*  from  trao  north. 

**  Monsvrod  from  tiw  tail. 
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<» 


I 

Radius  of  ths  Earth 

Mjatnitoni 

ft 

1 

! 

f 

S 

Rafaranca  araa  of  tha  misslla 

ft^ 

1 

( 

S« 

C  (1  -  f^)  (Ref  6) 

ft 

! 

T 

Thrust  vector  of  tha  missile 

i 

i 

f 

t 

Tima  -  indapandant  variable 

sec 

1 

1 

"^AS 

Rafaranca  temperature  of  tha 

f 

atmosphere 

J 

i 

Ta 

Temperature  of  the  atmosphere 

°K 

t 

i 

i 

’^Tt 

Thrust  known  for  input  data 

i 

Ttv 

Vacuum  thrust  of  the  missile 

l^f 

i 

f 

1 

V 

Velocity  of  the  missile 

ft/  sec 

f 

^SD 

Speed  of  sound 

ft/  sec 

o 

^SDS 

Reference  speed  of  sound 

ft/ sec 

^x 

Velocity  components  along  X,  Y,  Z 

ft/  sec 

1 

Vwx*V*V 

Wind  velocity  components  along 

X,  Y,  Z 

ft/ sec 

[ 

f 

J 

Vi.  V, 

Velocity  components  along  axis  1  k  2 

ft/ sec 

1 

j 

X,  Y,  2 

Geocentric  coordinate  system 

ft 

I 

1 

[ 

*L*  ^L* 

Launch  coordinate  system 

ft 

1 

i 

Xj.X^.Xj 

Missile  nonrotating  coordinate  system 

None 

[ 

Xj.X^.x'j 

Missile  rotating  coordinate  system 

None 

a 

GREEK  LETTERS 

Angle  between  east  and  axis 

deg 

Thrust  misalignment  angle 

rad 

i 

1 

e 

Euler  angle  aa  shown  in  Rgure  2 

deg 

®c 

Geocentric  latitude 

1 
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PimgBitani 

Geodetic  latitude 

deg 

Density  of  the  atmosphere 

slugs/ ft^ 

Summation  sign 

None 

Euler  angle  as  shown  in  figure  2 

deg 

Angle  at  which  the  thrust  misalignment 
acts  in  the  X2,  X^  plane 

rad 

Angle  at  whi^h  th^  thrust  misalignment 
acts  in  the  X2,  X^  plane 

rad 

Angular  velocity  of  missile  in  the 

X^,  X2,  X^  coordinate  system 

rad/ sec 

Angular  vejocity  of  the  missile  in  the 

Xj^,  X2,  Xj  coordinate  system 

rad/ sec 

s  u^k  =  Angular  velocity  of  the  Earth 

rad/ sec 

SUBSCRIPT  REFERS  TO 


C 

E 

G 

i 

L 

P 

S 

T 

V 

X 

X.  y.  » - 
w 

1,2,3 


a 


Geocentric 
The  Earth 
Geodetic 
Initial 

Launch  system 
Propellant 

Reference  conditions 
Thrust 
Vacuum 
-  Axis 

X,  Y,  Z  Coordinate  system 
Wind 

X^,  X2,  Xj  Coordinate  system 
Angle  ol  attack 
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An  arrow  over  a  variable  Indicates  that  the  variable  Is  a  vector. 

A  dot  over  a  variable  Indicates  the  time  derivative  of  that  variable. 


TDR-63.U 


(1)  Coordinate  eyeteme. 

(a)  Earth-centered  coordinate  system  (X,  Y,  Z). 

A  right-handed,  Earth- centered,  Cartesian  coordinate 
system  is  used.  The  X,  Y,  and  Z  axes,  shown  in  figure  1  are  oriented 
as  follows: 

The  X  axis  is  the  node  of  the  equatorial  plane  and  the 
plane  containing  the  Z  axis  and  Greenwich  Meridian. 

The  Y  axis  also  lies  in  the  equatorial  plane  and  is  at 
right  angles  to  the  X  and  Z  axes. 

The  Z  axis  lies  along  the  spin  axis  of  the  Earth  and  is 
at  right  angles  to  the  X  and  Y  axes. 

(b)  Launch  coordinate  system  (X^  ,  Yj^  ,  Z^^). 

A  right-hand  coordinate  system  with  the  X^  ,  Y^  plane 
tangent  to  an  oblate  Earth  at  the  launch  point  is  used.  The  X^  ,  Y^^  ,  Zj^ 

axes  are  oriented  as  follows: 

The  X|^  axis  is  in  the  direction  of  the  Initial  launch 

aximuth. 

The  Yj^  axis  is  counterclockwise  90°  from  the  X^  axis. 

The  Zj^  axis  is  positive  along  the  local  geodetic  vertical. 

(2)  Development. 

(a)  Newton  second  law  for  a  rotating  coordinate  system. 

The  fundamental  equation  of  motion  for  a  particle  moving 
in  a  rotating  coordinate  system  such  as  the  Earth  is 

Ujj.  X  (ug  X  R)  (1) 

(Reference  t) 


d." 


2F 

M 


z/dl  X 

\dt 


t 


TDR-63-ll 


TDR-63-11 


where  R  ie  the  vector  distance  from  the  origin  of  the  rotating  X-Y-Z 
coordinate  system  to  the  particle, 

(R  =  IX  i  JY  -i-  kZ)  d  R  are  the  velocity  and  acceleration, 

^  '  dt  dt* 

respectively,  of  the  particle  measured  with  respect  to  the  rotating  axes, 

♦  ♦  •  ♦  •  ♦  •  ^  •  * 

iS.  =  iX  +  jY  +  kZ  and  sl3  =  iX  +  jY  +  kZ.  (2) 

dt  ^ 

Ujg.  Is  the  angular  velocity  of  the  Earth  (axis  system)  and  is  along  the 
Z  axis,  (u  =  ku)|^ ).  2?  Is  the  vector  sum  of  the  forces  acting  on  the 

particle:  thrust,  drag,  and  gravity,  (2F  =  ^ 

the  total  instantaneous  mass  of  the  particle.  ^  coriolls 

pseudo-acceleration  of  the  axis  system  due  to  the  rotation  of  the  Earth; 

x  (u^  X  R)  is  the  centrifugal  pseudo-acceleration  of  the  axis  system  due 

to  the  rotation  of  the  Earth. 

I  t 

3  ^ 


llB  X  Up 

dt  ^ 


I 

dt 


iZ 

dt  dt 
0  1 


2(f  «  "El  .  I 


(u^  X  R) 


dt 


dt 


0 

X 


0 

Y 


1 

Z 


u)p(  -  lY  +  jX) 


“e  X  (u j,  X  R) 


1 

0 

i-Y 


0 

X 


1 

0 


up^(  -  IX  -  jY) 


(3) 
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From  oquattons  2  and  3  the  components  of  equation  1  can  be  expressed: 


d^X 

dt2 

S  2  Fx 

M 

dX  +  Ug^X 
dt 

dt^ 

=  j£Zy 
M 

- 

dX  + 
dt 

d!z 

=  ZFa 

dt^ 

M 

The  program  assumes  that  there  are  three  forces  acting  on  a  rocket 
vehicle:  thrust,  drag,  and  gravity.  Thrust  is  sasumed  to  act  along  the 
longitudinal  axis  of  the  rocket  vehicle.  Drag  is  assumed  to  act  opposite 
the  velocity  vector.  Gravity  is  assumed  to  be  directed  toward  the  center 
of  mass  of  m  oblate  Earth. 

(b)  Gravitational  attraction  equations. 

The  equations  of  gravitational  attraction  are  obtained 
from  reference  8,  and  converted  to  a  Cartesian  coordinate  system.  They 
are  as  follows: 

G  «  -GM  1  aK  -  15 

G  •  *  JB  -  15 


G.  -  -CM  ^—ll  +  JK.  -  15 

*  r"L  r’ 

Finally,  the  equations  for  total  acceleration  components  along  X,  Y,  and 
Z  axes  can  be  written  as  follows: 
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X  »  Is.  -  EL 

M  M 


Y  =  ly  -  C. 
M  M 


it 


+  2uj.Y  +  Uj.  X 


Gy  -  2«J.X  +  Uj.‘Y 


where 


Z  =  Is  -  H 

M  M 


-  G 


jc^  +  Y^  + 


(X 


ee  **9  **9  tA 

R  *  (X*  +  Y^  + 


(6) 


(7) 


(c)  Thrust  and  mass. 

The  mass  M  of  the  rocket  vehicle  and  the  thrust  T 


Tt 


at  air  pressure  are  assumed  to  be  known  functions  of  time.  In  general 

these  functions  will  be  nonlinear  and  discontinuous.  The  thrust  T  corres¬ 
ponding  to  air  pressure  P  nug  be  computed  from  the  equation 

ft 


T 


Tt 


+  (P 


at 


-  P  )A 

a  e 


(8) 


where  is  the  nozale  exit  area.  (A^  is  measured  in  square  inches; 
therefore,  P^^  and  P^  are  absolute  pressures  in  Ib/in^. ) 

In  this  trajectory- computation  program,  the  thrust 
function  is  evaluated  by  table  look-up  and  interpolation.  Since  this  function 
is  discontinuous,  a  set  of  tables  is  needed  for  each  stage  of  the  rocket 
vehicle.  The  tables  in  the  /orking  storage  are  changed  each  time  a  stage 
is  dropped. 

The  thrust  function  T.^^  vs.  t  is  given,  but  the  mass 
function  M  vs,  t  is  computed  by 
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where 


Tv 


Tt 


at 


A 

e 


(9) 

(10) 


Subscript  i  indicates  values  at  ignition,  b  denotes  values  at  burnout,  and 
Mp  is  the  mass  of  the  propellant.  For  clarity,  subscript  n  denoting  the 
corifiguriitioii  number  has  been  omitted  from  the  subscripted  symbols. 

(d)  Aerodynamic  drag. 


equation 


The  aerodynamic  drag  may  be  computed  from  the 
D  =  l/i  p  SV^Cj^(M.N.  )  (11) 


where 


p  =  Density  of  the  air  about  the  rocket  vehicle,  slugs/ 
ft^ 


Cp  (M.  N.  )  =  Drag  coefficient,  assumed  to  be  a  function 
of  Mach  number  only  in  this  report 
(dimensionless) 

S  =  Cross  -section  area  of  the  rocket  vehicle,  ft^ 

M.  N.  =  Mach  number  =  ^/^SD 

Vgp  =  Velocity  of  sound  in  the  air  about  the  rocket 
vehicle,  ft/ sec. 


It  should  be  noted  that  each  stage  of  the  rocket  vehicle 
could  have  a  different  diameter.  The  diameter  of  the  largest  stage  in  the 
assembly  that  has  not  dropped  off  will  be  used. 

Two  drag  coefficients  for  each  configuration  of  the 
rocket  vehicle  are  needed.  It  is  assumed  that  the  drag  of  a  configuration 
during  powered  flight  is  different  from  (less  than)  the  drag  before  ignition 
or  after  burnout  by  the  amount  of  the  base  drag.  The  drag  coefficients  for 
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the  powered  and  coasting  conditions  will  be  denoted  by  and  , 

respectively,  with  numerical  subscripts  added  to  denote  the  configuration  or 
stage  number.  Thus,  up  to  six  different  drag  coefficients  will  be  needed  for 
a  three- stage  rocket  vehicle:  f  '  *^DB3  ’  ^DCl  ’  ^DCZ  * 

^DC3  ■  ^DCI  s^>id/or  ai’c  not  needed,  of  course,  if  the  first  and/or 

second  configuration  do  not  coast. 

The  neglect  of  yaw  angle  in  determining  the  drag  co¬ 
efficient  is  justified  by  the  fact  that  the  thrust  will  dominie  the  motion  so 
that  fairly  large  errors  in  the  aerodynamics  will  have  little  effect  on  the 
trajectory.  Also,  the  yaw  angle  (and  its  effect  on  the  coefficients)  will 
presun'iafciy  be  small. 

(e)  Atmosphere  and  wind. 


The  values  of  p 
15 


SD  ' 


and  P  are  available  from 
a 


the  COESA  1962  model  atmosphere  or  from  launch  site  soundings. 
Actually  is  not  measured  directly;  instead,  the  air  temperature  T 

is  recorded  and  is  computed  from  the  equation 


SD 


(T  ,  T 

SDs  a /  as 


(12) 


where  is  the  standard  velocity  of  sound  (ft/ see)  corresponding  to 

standard  air  temperature  T  ,  and  the  units  of  T  and  T  should  be 

as  a  as 


If  a  wind  is  blowing,  it  will  have  an  important  effect  on 
the  trajectory  of  a  multistage  unguided  rocket  vehicle.  The  effect  is  most 
pronounced  during  the  first  stage,  and  decreases  thereafter.  After  burnout 
of  the  last  stage,  the  effect  will  be  fairly  small  and  can  be  neglected.  The 
wind  may  be  taken  into  account  by  computing  the  velocity  components  from 
the  following  equations; 


(13) 
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axis  systam.  Th«  total  velocity 
computed  from 


the  wind  components  In  the  Earth-centered 
V  with  respect  to  the  air  mess  may  be 


V 


•4  ♦ 


(14) 


Equations  13  follow  from  the  definition  of  ,  Vy  ,  and  as  com¬ 
ponents  of  the  velocity  of  the  rocket  vehicle  with  respect  to  the  air  mass. 


It  is  assumed  that  the  wind  velocity  is  a  function  of 

only  and  is  horisontal  (that  is,  ~  wind  velocity  is 

L 

exactly  horizontal,  then  Vy^2  f  ^  ^XY  ^  ^WZ  ~  ^ 

Li  Li  Li 


may  be  thought  of  as  a  flit-earth  approximation,  but  unless  very 

large,  it  is  doubtful  whether  can  be  measured  accurately  enough  for  the 

approximation  to  be  questionable. 


The  meteorological  data  taken  before  a  launch  should 

include  the  wind  velocity  and  the  wind  azimuth  angle  as  well  as 

p ,  ,  and  vs.  altitude.  By  convention  Ay^  is  defined  to  be  the 

azimuth  angle,  measured  clockwise  from  the  North,  from  which  the  wind 

is  blowing;  then  Ay^  is  also  the  azimuth  angle,  measured  clockwise  from 

the  South,  to  which  the  wind  Is  blowing.  It  will  be  seen  that  Vy^^  and 

L 

computed  from  the  following  equations: 

WY^ 


cos  (Ay^ 

-  Ay  ) 

fin  ) 

Lj 

(15) 

^WX  * 


and  V 
^WY 


y^y  are  then  rotated  into  the  Earth- centered  components 
and  * 

The  presentation  of  ttw  linear  accrteratlon  eqpMitloiie  for 
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the  burning  period  bee  now  been  completed. 


(1)  Coordinate  eyeteme. 

(a)  Nonrotating  body  axis  ,  X2  ,  X^) . 

A  right-hand  Cartesian  coordinate  system  which  is  fixed 
to,  but  does  not  spin  with  the  body.  The  X^  axis  is  along  the  spin  axis. 

The  X^  axis  lies  in  the  vertical  plane,  while  the  X^  axis  lies  in  the  hori¬ 
zontal  plane.  This  axis  system  is  related  to  the  reversed  launch 

coordinate  system  by  the  two  Euler  angles  6  and  q> . 

(b)  Rotating  body  axis  (X^,  X^,  X^). 

A  body  fixed  coordinate  system  with  the  same  origin  as 
the  Xj^,  X2.  X^  system.  Axes  X^  and  X^  are  along  the  transverse 
principal  axis  of  the  rocket  vehicle. 

(2)  Development. 

(a)  Body  axis  equations. 

The  angular  acceleration  equations  of  the  rocket  vehicle  may 
be  developed  from  the  vector  equation  of  reference  1.  The  only  torque  con¬ 
sidered  is  the  overturning  moment.  All  other  moments  including  pitch  and 
jet  damping  are  assumed  to  be  smdl  and  will  be  neglected;  the  rocket  vehicle  is 
spun  to  reduce  the  effects  of  any  misalignments,  asymmetries,  or  unbalances, 
but  such  a  spin  is  not  large  enough  to  develop  an  appreciable  Magnus  nooment. 
The  vector  equation  of  angular  motion  is  then  as  follows: 


dS  =  dH 
dt  dt 


0  X  H  X 


where 


Angular  nunnentom  vector,  Ib-ft-sec 
Overturning  moment  (vector),  Ib-ft 

Angular  volocity  of  the  X^,  X2,  system,  rad/ sec 
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Th*  angular  valocity  of  tha  miaaiia  ia  doaotad  by  u  ;  tKa 
componanta  of  u  in  tha  X^,  ayatam  ara  -t-  NtQ2»fi3.  whare  N 
ia  tha  axial  apin  of  tha  miaaiia  and  0^^,  Qj  ara  tha  componanta  of  ?  . 
Than  aquationa  for  "S  and  H  ara  aa  followa:  (rafaranee  1. ) 


whara 


Q  a  IfOl  ♦  l^Oj 

H  aTj  A(Qj  +  N)  +  T2  BQ2  +  Tj  BQj 

(17) 

A  a  noomant  of  inartia  of  the  rocket  vaUcla  about  tha 
longitudinal  principal  axla,  alug>ft^ 

B  a  moment  of  inartia  of  tha  rocket  vehicle  about  a  , 
tranavaraa  axia  through  tha  center  of  maaa,  alug-ft^ 


It  ia  aaaumed  that  tha  maaa  diatribution  ia  aymmetric,  ao 
tha  moment  of  inertia  ia  the  aame  about  any  tranaverae  axia  through  the  canter 
of  maaa. 


(b)  Momenta  of  inertia. 

It  ia  aaaumed  that  an  intemal-buming  aolid  propellant  ia 
uaed.  Approximate  valuea  of  A  and  B  during  burning  are  computed  by  uae 
of  the  following  formulae: 


A  = 

^  - 

kip 

(18) 

B  r: 

®i  ' 

Mj  (G  -  Gj)  (Gj  -  Gp)  -  (Mj  -  M)kjp 

(19) 

whare 

G  = 

(Gj  -  Gp)  (Mj  -  M)/M 

(20) 

Ae  before,  aubacript  i  indicataa  valuae  at  time  of 
ignition  and  P  danotaa  a  property  of  tha  propellant;  i*  the  radiua  of 

gyration  of  tha  propellant  grain  about  ita  longitudinal  axia  (ft),  Kgp  !•  iha 
radiua  of  gyration  of  tha  propallaitt  grain  about  a  tranavaraa  axia  through  ita 
canter  of  maaa  (ft),  and  G  ia  tha  diatanca  from  tha  baaa  to  tha  caittar  of 
maaa. of  the  rocket  vehicla  (ft).  Tha  quantitiaa  Gp  ,  k^p  and  kg p  ara 
aaaumed  to  be  conatant  for  an  internal  burning  grain;  A^,  B^.  G^,  and  M^ 
ara,  of  couraa,  conatam,  ao  M  ia  tha  only  variable. 
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Eqaationf  19  and  20  cannot  be  naed  for  an  end-burning 
grain  unless  formulas  are  added  for  and  Gp  .  None  of  these  equations 

apply  for  liquid  propellant  rockets. 

(c)  Aerodynamic  moment. 

The  expression  for  is  as  follows: 

-  V4pdSVC„,  (TjVj  -  TjVj)  (21) 


For  the  present  application  is  the  static  stability 

derivative  (dimensionless)  and  is  a  negative  number  which  can  be  computed  by 
use  of  the  equation: 


'Ma 


(22) 


s  Normal  force  coefficient  (dimensionless) 

3  Distance  from  the  base  of  the  rocket  to  the 
normal  force  center  of  pressure,  feet 

In  this  report  ,  and  ,  are  assumed  to 

be  functions  of  Mach  number  only.  In  the  trajectory- computation  program 
these  functions  and  are  evaluated  by  table  look-up  and  interpolation. 

(d)  Euler  angle  relation. 

If  equations  17  and  21  are  substituted  into  equation  16, 
the  and  components  of  the  resulting  vector  equation  will  be  as  follows: 

=  Msp  VdS 

BQj  +  BDj  +  (B  -  A)  -  ANQ^  =  -V^pVdS  (23) 

It  is  desirable  to  replace  in  these  equations  by 

functions  of  f  and  6  .  This  may  be  done  by  use  of  the  following  relations: 

Qj  3  -  f  sin  6 


where 


'Na 


18 


The*e  equations  neglect  the  angular  velocity  of  the  Earth  (u^);  they  are 
written  by  inspection  of  figure  2.  It  will  be  seen  that 


Figure  2.  Coordinate  system  for  angular  equations. 
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Substitution  of  equations  24  and  25  into  equations  23  yields  the  required 
angular  acceleration  equations  as  follows: 

-B0  -  B0  +  (A-B)(p^  sin  0  cos  G  -  AN  qi  cos  0 
=  'ApVSdCj^^  V3 


-Bq)  cos  0  -  Bip  cos  0  +  (2B  -  A)  q)  0  sin  0  +  AN  0 


=  -ttpVSdC„,  Vj 


(261 


These  equations  are  dynamically  exact,  except  for  the 
neglect  of  which  is  negligibly  small  compared  to  q>,  0,  and  N. 

Equations  26  may  be  put  in  a  form  more  suitable  for 
computation  by  dividing  through  by  B,  The  equations  become 


*0  +  (B/B)  0  +  (1 
=  (-VipVdSC^^) 


-  (A/B)  )  q)^  sin  0  cos  0  +  (A/B)N  ^  cos  0 


q)  cos  0  +  (B/B)  9  cos  0  -  (2  -  (A/B)  )  q)  0  sin  0  -  (A/B)N0 
VV, 

=  Wp  B 

(27) 

(e)  B  Terms. 

For  an  internal  burning  solid- propellant  rocket,  the 
following  formulas  are  used  to  compute  the  B  terms: 


B  =  M 


where 


These  formulae  are  easily  derived  from  equations  9  smd  19. 


(28) 
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(f)  Body-crot«  wind  components. 

A  coordinate  transformation  is  needed  to  compute 

and  Vj  from  ,  Vy  ,  .  The  following  equations  can  be 

jLi  Li  Li 

written  from  inspection  of  figure  2: 


-  V„  cos  9  + 
*L 


V9  sin  cp 


=  -Vv  eln  7  sin  0  -  Vv  cos  6  -  V_  cos  ^sln  0 

(29) 


This  completes  the  derivation  of  the  equations  of  motion  during  burning, 
c.  Summary  of  cqttatiQBi  of  motion. 


If  the  rocket  is  assumed  to  have  thrust  misalignments,  additional 
terms  are  added  to  the  equations  to  account  for  the  new  forces  and  moments 
created.  For  convenience,  the  equations  are  restated  here  with  the  new 
terms  added. 


X  = 


Y  = 


Z  = 


i 

M 


1_ 

M 


M 


+  G 


+  G 


2u>,-  Y  + 

u,2 

X 

E 

£ 

2uj.X  + 

4 

Y 

tx  -  4] 
k  -4] 

k  -  4] 

"  ^  9  ♦  {[(2  -  ^  )  9  sin  0  +  ^  nJ  0 

1C„-G)  TCtTCO.»A1  _L_ 

8  Na  P  B  B  /  cos  0 

"  B  ®  ”  1^1  -  9  Bln  8  +  ^  p  cos  0 


VV-  _  T  G  fc  ™  sin  9 


(30) 
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where 

&  thrust  misalignment  angle,  radians 

q>rp=  orientation  angle  of  jet  misalignment  force;  measurea  in  the 

-  X^  plane  from  the  X^  axis  and  positive  in  the  sense 

of  a  -.ounter  clockwi£>e  rotation  as  seen  from  the  positive  X^ 

axis,  ra<'ians. 

t 

9^  -p  +  J  measured  in  the  X^  -  X^  plane  from  the 

ti 

X^  axis  and  positive  in  the  sense  of  a  counterclockwise 
rotation  as  seen  from  the  positive  X|  axis,  radians. 

These  equations  require  a  table  of  roll  rate  (N  =  vs.  *.ime. 

The  thrust  vector  is  given  in  the  launch  coordinate  system  and 
rotated  to  the  Earth- centered  system  by  the  use  of  a  matrix  rotation.  The 
thrust  vector  in  the  launch  coordinate  system  is 


T^j^  =  T  j^cos  9  sin  9  -  b.p(cos  9^^  cos  9  + 
sin  9^  sin  9  sin  0  )j 

Tyj^  =  -T  ^sin  6  +  6  .j,  sin  9^  cos  oj 

T^^  =  T  ^cos  9  cos  9  +  6  .j,  (cos  9^  sin  9  -  sin  9^  cos  9  sin  0) 

(31) 

d.  Geocentric  relationships. 

(1)  Shape  of  the  Earth. 

Figure  3  shows  a  meridian  section  of  the  earth  where  a^  is 
the  semi-major  (equatorial)  axis,  9^  is  the  geocentric  latitude,  and  is 
a  radius  vector  from  center  to  the  surface  of  the  earth.  R^  is  a  function  of 
the  geocentric  latitude  and  is  given  as 
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Hq  -  Geodetic 
Altitude 
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To  convert  f  odotlc  lotttudo  and  altitude  to  geocentric  latitude  and  radius: 
(reference  6) 


tan  6 

c 


LL^ 


tan 


(33) 


where 


C  A 


iE. 


(I  -  4  .1."  6^)'^ 


S  4  C(1  -  e“) 


R  = 


[(C*  H, 


coe^  +  (S  +  Hq)^  eln^  0^ 


yi 


(34) 


(35) 


to  convert  geocentric  latitude  and  radius  to  geodetic  latitude  and  altitude: 
(reference  7) 


Gq  =  e 


»G  =  ■  *E 


[t  »in  2  ©c  +  4  ec(-^  - 

j^-f  .in=^  Gc  2  -4^ 


(36) 

The  following  geocentric  constants  were  used  in  the  program. 
Adopted  Geocentric  Constants  (1961)  (reference  8) 


aj,  =  20,  925.  647.  IE  ft 

GM  =  1.4076427  X  10^^  ftV»ec^ 

=  2.  316686  X  10^^  NM  (^  )^ 

sec 

gj,  =  32.  174  ft/sec^ 

1/f  a  298. 30  t  0. 05 
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O 

a  0.0818133302 
J  =  {1623.42  ±  0.5)  X  10’^ 

Pj.  =  0.00673852 

Ug  =  7.292115  X  10"®  rad/ sec 
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INPUT  DATA 

All  input  data  are  read  into  SPURT  at  one  time  and  in  one  "read"  block. 
This  includes  data  for  the  subroutines  (two-body  and  impact).  Most  data, 
except  for  special  cases  and  control  integers,  are  read  in  a  FORTRAN 
7F10.  0  format.  This  format  is  for  seven  data  words  (of  ten  digits  each)  per 
card.  The  decimal  point  is  always  punched  in  the  field.  On  Repeat  Blocks, 
as  many  words  will  be  read  as  specified  by  a  previously  read  integer. 
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nan  (mm  or  imis  im: 


PazaiBBtcr 

Dimensian 

Cbluan 

Mode 

Bemsrks 

No.  of  stages 

none 

1 

I 

■■■ 

1  Card 

Payload  velt^ 

pounds 

KBS 

V 

1 

1  Saae 

1-80 

Hollerith 

1  Card  1 

Initial  latitude 

degrees  (tN) 

1-10 

V 

i 

1 

Initial  longitude 

degrees  (t  E) 

11-20 

V 

Initial  altituie 

feet 

21-30 

V 

degrees  CN  from  N 

31-40 

V 

1  Card 

Initial  X 

feet 

1 

o 

V 

Initial  Y 

feet 

51-60 

V 

Initial  Z 

feet 

6l  -  70 

V 

Initial  Si 

feet/sec 

1-10 

V 

k 

Initial  Y 

feet/see 

11-20 

V 

Initial  Z 

feet/sec 

21-30 

V 

Initial  9 

degrees 

31-40 

V 

1  0 

ard 

Initial  e 

degrees 

o 

UN 

V 

Initial  7 

deg/sec 

51-60 

V 

Initial  e 

deg/seo 

61-70 

V 

Start  tlae  >  0. 

seeoDds 

1 

H 

V 

1  C 

1 

ird 

Time  of  last  stage  B»0. 

seooBds 

11-20 

V 

_j 

' 

© 
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lASter  co:atrol  l 


Control  1 

none 

1  -  2  . 

I 

metro  &  wind 

Control  2 

none 

3  -  4 

I 

print  Input  data 

Control  3 

none 

5-6 

X 

vrite  BAIT  tape 

Control 

none 

7  -  3 

I 

two  'body  sub. 

Control  5 

none 

9-10 

I 

linea/page 

Control  6 

none 

11-12 

T 

•tm 

’.nrite  plot  tape 

Control  7 

none 

13  -  14 

1 

rl^t  hand  co3:d. 

Control  3 

none 

15  -  l6 

I 

inrpact  sub. 

Control  9 

17  -  18 

I 

“all  thl^ 

Control  10 

none 

19  -  20 

I 

"blccl:  on 

iro.  of  sn;'.n  nal'.'.Gs 

none 

21  -  22 

I 

one  card 

Gv.tnv.t  t;r"0.  I'o. 

none 

23  -  2'-r 

I 

Do'j.  0*  intsr-.'o.lo.tion 

none 

25  -  2G 

I 

3_,C:’  tC..i'3  tr.’:le 

[  cecends 

7010.0 

[  11  spin  Yolues 

1 

i*cd/oec 

1  spin  TBlues 

Hi 

Gta'jc  nv:l'oor 

none 

1-2 

I 

iro 

Size  0?  thrust  table 

none 

3  -  4 

I 

N  (1,J) 

Size  of  Cq-  table 

none 

5  -  6 

I 

H  (2,J) 

Size  of  Cq^i  table 

none 

7-8 

I 

N  (3,J) 

Size  of  Cp  table 

none 

9-10 

I 

N  (4,J) 

Size  of  C|Jq  table 

none 

11-12 

I 

N  (5#J) 

© 
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Pressure  at  tlirust  meas. 

o 

pounds/in 

1-10 

Exit  area 

inch^ 

11-20 

Staco  cLiaiueter 

feet 

21-30 

lissile  C.G. 

feet 

31-40 

Dta~e  fuel  C.G. 

feet 

41  -  50 

Euc.l  a;:ial  I/'E 

feet^ 

51-60 

.\icl  trruiGTCrse  I/M 

feet^ 

61  -  70 

iissile  a;d.al  I 

olug-ft^ 

1-10 

iiLssile  trunsverse  I 

slug-ft^ 

11-20 

StaGC  ;;eiGht 

pounds 

21-30 

otaje  consuaed  \rt 

pounds 

31-40 

1 

Orientation  of  I.M.A. 

radians 

41  -  50 

Thrust  aisalicn  angle 

radians 

51-60 

-  Oi-ET  - 

61-70 

Ignition  time  from  lavmch 

seconds 

1-10 

Dumout  time  from  launch 

seconds 

11-20 

TICC  from  launch 

seconds 

21-31 

1  Card 


1  Card 


1  Card 


Repeat  fran  A' 
for  each  stage  (RS) 


m 
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WPP  AMP  maKtaum  npw?  oa 


Skip  Metro  Data 


No.  of  altitude  'Values 


Altl-tud-e  'values 


d 


' 

YES 

None 

1-2 

1 

KMEIBO 

feet 

7I10.C 

V 

WAIX 

' 

— <^QrERL  ] 

' 

HO 

r 

°c 

7710.0 

mm 

No.*  lOffllBO 

Kcai/H^ 

7710.0 

D 

Milliljars 

MSSSMMMMSaiEMM 

KMEERO 


feet/sec 


7P10.0 


No.* 


T11I-63-IL 


Tim-63-11 


1  ~  IMEBOSt  -  Hxitd  point  nunbors 
V  -  ^ASUMR  •  Floating  point  »ag>arB 
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MASTER  COMTRQL  NUMBEBS 

These  numbers  control  various  options  of  the  program  as  follows: 

Control  1.  When  KNTRL(l)  =  1*,  the  program  will  read  in  the 
temperature,  pressure,  and  density  and  use  them  to  compute  the  aerodynamic 
forces  and  moments.  Otherwise,  the  standard  1962  atmosphere  is  used. 

When  KNTRL(l)  =  1  or  2,  the  program  will  read  in  wind 
velocity  and  wind  azimuth  vs.  altitude.  The  program  will  compute  the  path  of 
the  missile  due  to  the  winds. 

Control  2.  When  KNTRL(2)  =  1,  the  input  data  are  printed  out. 

Control  3.  When  KNTRL(3)  =  N,  the  program  will  write  tape 

number  N  for  input  to  the  BATT  program. 

Control  4.  When  KNTRL(4)  *  1,  the  program  will  use  the  Keplerian 
(TWO-BOD)  subroutine  for  coasting  flight  on  the  last  stage.  If  KNTRL(4)  =  Z, 
it  will  use  it  on  all  stages. 

Control  5.  When  KNTRL(5)  =  NO,  the  program  will  write  NO  lines 
of  output  per  page. 

Control  6.  When  KNTRL(6)  =  N,  the  program  will  prepare  an  input 
tape  (N)  for  the  plot  program. 

Control  7.  When  KNTRL(7)  =  1,  the  launch  coordinate  system 

printed  out  is  a  right-hand  one  instead  of  a  left-hand  one. 

Control  8.  If  KNTRL(8)  =  1,  the  integrating  unpowered  flight 
trajectory  is  computed  for  all  stages  except  the  last  one.  If  KNTRLfS)  s  2, 
the  trajectory  is  computed  for  all  stages. 

Control  9.  Not  used, 

CantroLU.  Not  used. 


*  When  a  KNTRL  number  is  left  blank,  that  option  will  be  ignored. 
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INPUTS 


AALT 

AAZIM 

AAZI2 

ACODE 

ACODES 


ALAT 

ALA2 

ALON 

AL02 


AXLM 

AXLMB 


BAZIM 

BOOT 


BDOTB 


BLON 


CALAT 

CBOOCK 

COM 


Initial  axial  moment  of  inertia  for  each  etage  (slag.ft^) 
Initial  geodetic  altitude  (ft) 

Initial  launch  aaimuth  (deg) 

Initial  launch  aaimuth  (rad) 

Integration  code 

Integration  code  table 

Exit  area  of  rocket  for  each  stage  (in^) 

Initial  geodetic  latitude  (deg) 

Initial  geodetic  latitude  (rad) 

Initial  longitude  (deg) 

Initial  longitude  (rad)  4- 

Geocentric  altitude  (ft) 

2  4-  longitude  (rad) 

Colatitttde  (rad) 

Axial  moment  of  inertia  (elug.ft^) 

Axial  moment  of  inertia  over  traneverse  moment 
of  inertia 

Initial  transverse  moment  of  inertia  for  each  stage 
( slug-ft^) 

AAZ12.^  (rad) 

Rate  of  cjiange  of  the  transverae  moment  of 
inertia,  B  ( slug. ft2/ sec) 

B/  B  (sec"^) 

Integration  data  storage 

^  4-  AL02(rad) 

i  REE/^1  .  e^  Sia^  0^  (ft) 

Cosim  of  the  initisd  latitude 

CommM  block  for  integration 

CjyA/m  .  drag  parameter  used  lor  empty  stages 

(ft^/clttga)  * 


*  gMiia4'6lfb 

'J; 
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CHECK 

Flag  for  integration  check 

CN  * 

Input  table  of  normal  force  coefficient  C^  (rad**^ 

(X 

CNI 

Normal  force  coefficient  used  in  calculation 

CNMACH  * 

Input  CN  Mach  table 

CNMTAB 

Inverted  CN  Mach  table 

CNTAB 

Inverted  CN  table 

COEF 

Coefficient  used  in  aerodynamic  moment  = 

■> 

(Cp.CG)J^^ 


a 


COLAT 

Initial  co> latitude  =  '2  ~  ALA2  (rad) 

COTHC 

Co- geocentric  latitude  =  ^  ~  THC  (rad) 

CP 

* 

Initial  center  of  pressure  table  (feet  from  tail) 

CPHl 

Cosine  of  PHI 

CPHT 

Cosine  of  PHT 

CPI 

Center  of  pressure  used  in  calculation 

CPMACH 

* 

Initial  Mach  No.  table  for  center  of  pressure 

CPMTAB 

Inverted  CP  Mach  No.  table 

CPTAB 

Inverted  CP  table  (ft) 

CTHC 

Cosine  of  initial  geocentric  latitude  =  Cos  (THC) 

CTHET 

Cosine  of  THETA 

D 

* 

Diameter  of  each  stage  (ft) 

DIMACH 

Input  Mach  No.  table  for  burning  drag  coefficient 

D2MACH 

* 

Input  Mach  No.  table  for  coasting  drag  coefficient 

DC 

Direction  cosines  matrix 

DELPR 

* 

Print  times  used  in  impact  (sec)  + 

DELTH 

Difference  between  6  and  6  (rad) 

c  g 

DM  TAB 

Inverted  drag  Mach  table 

DRAGl 

* 

Input  table  of  burning  drag  coefficient 

DRAG2 

* 

Input  table  of  coastii^  drag  coefficient 

DT 

Print  interval  (sec) 

DTAB 

Inverted  drim  ts^le 

-f  Stored  ia  cooHOMm 
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1  • 

iMEttlS 

i  „  - . 

i 

DUM 

* 

Dummy  variable 

f 

t 

ERR 

Integration  error 

FDRAG 

Drag  force  (lb) 

FORCE 

Vacuum  thrust  (lb) 

FT 

Thrust  on  vehicle  (Ib) 

FTT 

Total  impulse  (lb  -  sec) 

G1 

Center  of  gravity  at  any  given  time  (feet  from  tail) 

GK 

Earth  oblateness  term  a^  J  (ft^) 

GM 

2 

Gravitational  constant  for  the  Earth  (ft  /sec  ) 

GO 

* 

Initial  CG  of  the  remaining  missile  for  a  given  state  (ft) 

GOP 

Difference  betv7een  GO  and  GP  (ft) 

GP 

* 

Initial  CG  of  the  fuel  for  a  given  stage  (ft) 

GRAVO 

Gravity  constant  (32. 174  ft/sec^) 

GRV 

Term  used  to  compute  gravity  components 

GRVX 

Gravity  component  along  the  X  axis  (ft/sec^) 

GRVY 

Gravity  component  along  the  Y  axis  (ft/sec^) 

GRVZ 

Gravity  component  along  the  Z  axis  (ft/sec^) 

H 

Velocity  vector  (no  wind)  (ft/ sec) 

HALFPI 

X 

2 

I 

Utility  index 

ICODE 

Code  used  to  determine  integration  method 

i 

> 

U 

Utility  index 

III 

Utility  index 

INTN 

* 

Order  of  interpolation 

i 

IRA 

Utility  index 

; 

IT 

Print  table  count 

r 

J 

Utility  index 

1 

JP 

Lines  per  page  count  for  output 

t 

f 

j 

JTBl 

* 

Two  Body  input  (number  of  time  changes) 

+ 

1 

JTB5 

* 

Two  Body  input  (number  of  look>angle  stations) 

+ 

K 

Utility  index 

o 

KB  ATT 

Tape  number  for  BATT  tape  =  KNTRL(3) 

+ 

KIX 

Number  of  stages  for  subroutine  impact 

•¥  Stored  in  common 

1 
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KMETRO 

INPUTS 

* 

Size  of  Metro  tables 

KNTRL(IO) 

Controls 

KPLOT 

Tape  number  of  PLOT  tape  =  KNTRL(6) 

L 

Stage  count 

LSKIP 

Skip  aeiodynamics 

N 

Table  size 

NO 

» 

Input  card  check 

N1 

Thrust  table  size  for  different  stages 

N2 

Burning  drag  table  size  for  different  stages 

N3 

Coasting  drag  table  size  for  different  stages 

N4 

CP  table  size  for  different  stages 

N5 

CN  table  size  for  different  stages 

NAME 

* 

Page  title  (up  to  80  Hollerith  characters)  + 

NOE 

Number  of  equations  used  in  integration 

NOT 

* 

Output  tape  number  + 

NNX 

* 

Table  of  NRAT  values  + 

NRAT 

* 

Table  size  for  each  stage  in  impact 

NS 

* 

Number  of  stages  (NS  <  10) 

NSPIN 

* 

Spin  table  size  (NSPIN  <  100) 

NXIM 

♦ 

Number  of  stages  read  in  impact  drag  table 

OMEG 

Spin  rate  of  Earth  =  (rad/ sec) 

OMEG2 

2  2  2 

Earth  spin  rate  squared  =  (rad  /sec  ) 

OUT 

Variable  used  for  output 

PI 

Pressure  times  exit  area  (lb) 

PAT 

Pressure  at  which  thrust  is  measured  (Ib/in'') 

PHI 

Euler  angle  used  in  equations  (rad) 

PHID 

First  derivative  of  PHI  (rad/ sec) 

PHIDD 

Second  derivative  of  PHI  (rad/ sec  ) 

PHT 

♦ 

Orientation  angle  of  thrust  misalignment  (rad) 

PI 

It 

PRES 

Atmospheric  pressure  (Ib/ft^) 

PRTIME 

Time  at  which  to  print  (sec) 

+  Stored  in  common 
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PWGT  * 

PWGTC  * 

PX 

PXD 

PXDD 

PXL 

PXLD 

PXLDO 

rXND 

PYL  WGT  * 

R 

RA  ♦ 

RANGE 
RANGE 1 
RAD 

RB  ♦ 

RE 

REE 

REL 

RHO 

ROT 

ROTl 

RIX 

R2X 

R3X 

R4X 

R5X 

R2 


Stag*  input  weight  (lb) 

Stag*  fuel  weight  (lb) 

Earth  centered  position  vector  used  in  integration  (ft) 

First  derivative  of  position  vector  (ft/ sec) 

Second  derivative  of  position  vector  (ft/sec^) 

Position  vector  in  launch  coordinate  system  (ft) 

Velocity  vector  in  launch  coordinate  system  (ft/ sec) 

Acceleration  vector  in  launch  coordinate  system 
(ft/sec^) 

Velocity  vector  in  local  coordinate  system  (ft/ sec) 
Payload  weight  (lb) 

Distance  from  Earth  center  to  vehicle  (ft) 

Fuel  axial  radius  of  gyration  squared  for  each 
stage  (ft2) 

Range  at  burnout  of  each  stage  (N.M. )  + 

Dtimmy  variable  used  for  impact  (N.M . ) 

180 

Fuel  transverse  radius  of  gyration  squared  for 
each  stage  (ft2) 

Radius  of  Earth  as  a  function  of  latitude  (ft) 
Equatorial  radius  of  the  Earth  (ft) 

Distance  from  Earth  center  to  launch  pt  (ft)  + 
Atmospheric  density  (slugs/ft  ) 

Matrix  from  launch  to  Earth  centered  coordinates 
"CCMMON"  rotation  matrix  + 

Dummy  variable  used  in  X  integration 
Dummy  variable  used  in  Y  integration 
Dummy  variable  used  in  Z  integration 
Dummy  variable  used  in  9  integration 
Dummy  variable  used  in  6  integration 
Distance  squared  frtmi  Earth  center  to  vehicle  (ft^) 


S  4  C(1  -*^)(ft) 

SILAT  Sins  of  the  initial  latitude 

*  Stored  in  common 
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INPUTS 


SMAX 

Maximum  Integration  step  size  (sec) 

SMIN 

Minimum  integration  step  size  (sec) 

SPHI 

Sine  of  PHI 

SPHl 

Sine  of  PHT 

SPl 

Inverted  spin  table 

SPIN 

* 

Input  spin  table  (rad/ sec) 

SPIT 

Integrated  spin  table  ( rad) 

SPT 

Inverted  spin  time  table  ^ 

SPTIME 

* 

Input  spin  time  table  (sec) 

SS 

Integration  step  size  (sec) 

STARTT 

* 

Start  time  or  launch  time  (sec) 

STHC 

Sine  of  the  initial  geocentric  latitude 

+ 

STHET 

Sine  of  THETA 

STP 

* 

Initial  PHI  (deg) 

STPD 

* 

Initial  PHI  DOT  (deg/ sec) 

STT 

* 

Initial  THETA  (deg) 

STTD 

* 

Initial  THETA  DOT  (deg/ sec) 

STX 

» 

Initial  position  vector  in  launch  coordinate  system  (ft) 

STXD 

* 

Initial  velocity  vector  in  launch  coordinate  system 
(ft/ sec) 

SW 

Angle  between  wind  azimuth  and  X  axis  ( rad) 

T 

Thrust  vector  in  launch  coordinate  system  (lb) 

TABLE 

Table  of  printout  times 

TB2 

* 

Number  of  orbits 

+ 

TB3 

* 

Time  increment  At  (sec) 

+ 

TB4 

* 

Ending  time  (sec) 

+ 

TB6 

* 

Name  of  look- angle  station  (up  to  24 
characters) 

Hollerith 

+ 

TB7 

* 

Latitude  of  look- angle  station  (deg) 

■J' 

TBS 

* 

Longitude  of  look- angle  station  (deg) 

+ 

TB9 

* 

Altitude  of  look-  angle  station  ( ft) 

+ 

TDEL 

* 

Thrust  misalignment  angle  ( rad) 

TFIMP 

Stage  drop  time  (sec) 

+ 

THC 

Initial  geocentric  latitude  ( rad) 

■¥  Stored  in  common 
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INPUTS 


THETA 
THETD 
THE TDD 
THRUST  * 

THTAB 
TIME 

TIMET  * 

TMBO  * 

TMCC  * 

TMI  ♦ 

TOMEG 

TP 

TPl 

TP  2 

TPHO 

TRERR 

TRVM 

TSTOP  * 

TT 

TTTAB 

TWOPI 

V 

VA 

V2 

V3 

V  MACH 
VWX 
VWY 
VWZ 
VX 

WALT  * 


Euler  angle  in  the  horizontal  plane  (rad) 

First  derivative  of  THETA  (rad/ sec) 

Second  derivative  of  THETA  (rad/sec^) 

Input  thrust  table  (lb) 

Inverted  thrust  table 

Time  (dependent  variable)  (sec) 

Time  table  for  thrust  (sec) 

Burnout  time  for  each  stage  (sec) 

Time  to  change  coefficients  for  each  stage  <sec) 

Ignition  time  for  each  stage  (sec) 

Twice  the  Earth  spin  rate  =  2  Ug,  (rad/ sec) 

Intermediate  variables  used  in  output  block 

Intennediate  variables  used  in  output  block 

Intermediate  variables  used  in  output  block 

Orientation  angle  of  thrust  misalignment  at  a  given 
time  (rad) 

Integration  truncation  error 

2 

Transverse  moment  of  inertia  at  any  time  (ft  -  slug) 
Preset  time  to  stop  powered  portion  of  program  (sec) 
Thrust  vector  in  Earth  centered  coordinate  system  (lb) 
Inverted  thrust  time  table 

2it 

Velocity  of  vehicle  (ft/ sec) 

Velocity  of  sound  (ft/ sec) 

Cross  wind  perpendicular  to  the  velocity  vector 
(ft/ sec) 

Cross  wind  perpendicular  to  the  velocity  vector 
(ft/ sec) 

Mach  number  of  vehicle 

Wind  components  in  Earth- centered  system  (ft/ sec) 

Wind  components  in  Earth- centered  system  (ft/ sec) 

Wind  components  in  Earth- centered  system  (ft/ sec) 

Velocity  vector  with  wind  in  Earth-centered  system 
(ft/ sec) 

Input  metro  altitude  table  (ft) 
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USEllIS. 


WALTU 

Inverted  metro  altitude  table  (ft) 

WDEN 

* 

Input  metro  atmosphere  density  (kgm/m^) 

WDU 

Inverted  metro  atmosphere  density  (slugs/ft^) 

WGT 

Mass  of  remaining  missile  for  each  stage  (slugs) 

WGTAB 

Inverted  mass  table  (slugs) 

WGTC 

Mass  of  fuel  for  each  stage  (slugs) 

WGTCI 

Mass  expelled  from  ignition  (slugs) 

WGTI 

Mass  at  a  given  time  WGTI  =  WGT. WGTCI  (slugs) 

WlNDA 

* 

Local  wind  animuth  (deg  from  North,  clockwise) 

WINDV 

* 

Local  wind  velocity  (ft/ sec) 

WMAX 

Maximum  altitude  of  metro  table  (ft) 

WPRES 

* 

Input  metro  pressure  table  (millibars) 

WPU 

Inverted  metro  pressure  table  (Ib/ft  ) 

WTEMP 

m 

Input  metro  temperature  table  (**0 

WTU 

Inverted  speed  of  sound  table  (ft/ sec) 

wx 

Wind  component  parallel  to  Earth. centered  X  axis 
(ft/ sec) 

wy 

Wind  component  parallel  to  Earth.centered  Y  axis 
(ft/ sec) 

wz 

Wind  component  parallel  to  Earth.centered  Z  axis 
(ft/ sec) 

XIX 

Variables  used  in  geodetic  to  geocentric  conversion 

XIY 

Variables  used  in  geodetic  to  geocentric  conversion 

XDIMP 

Velocity  vector  at  burnout  for  impact  (ft/ sec)  -t- 

XIMPA 

Position  vector  at  burnout  for  impact  (ft)  + 

XJ 

First  harmonic  coefficient  for  oblateness 

XMACH 

Mach  number  table  for  Cjyk/m  + 

xsa? 

Spin  rate  at  any  given  time  (rad/ sec) 

XSPT 

Integrated  spin  at  any  given  time  ( rad) 

xvx 

Variable  used  in  calculation  of  PX 

Z2 

Lines.per.page  count  for  Two  Body  ou^ut  -f 

■f  Stored  1b  common 
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5. 


SPURT  FLOW  DIAGRAMS. 
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pure  SUM  o^sTAec 
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lOtl  I  YCS 


» 
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6.  SUBROUTINES. 
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a.  Atmosphere  subroutine 


The  1962  COESA  standard  atmosphere  is  incorporated  as  a  subroutine 
of  SPURT.  This  atmosphere  uses  the  equations  derived  for  the  1959  ARDC 
model  atmosphere.  These  equations  are 

(1)  Molecular- scale  temperature 


(2) 


Tm  (Tm)^  +  Lm  (H  - 
Pressure  -  altitude 


GM 


p=  p^ 


(TmK 


(Tm)^  +  Lm  (H  - 


R*Lm  for  Lm  j-  0 


P  =  P  exp 
b 


-GM^  (H  -  H^) 


R*  (Tm), 


for  Lm  =  0 


(3)  Density  -  altitude 


p  =  3. 236598  x 10 


■4  _JE_ 
Tm 


(4)  Speed  of  sound 

Vs  =  1116.4437  yfmTTnT 


where  b  =  subscript  that  refers  to  the  quantity  at  the  base  of  the  constant- 
gradiant  layer. 


GMo/R*  = 

constant  of  the  air  gas 

H 

geopotential  altitude  in  meters 

»b 

geopotential  altitude  in  meters  at  the  base  of  a 
constant  gradiant  layer 

Lm  = 

dTm/dH  =  the  gradiant  of  the  molecular  scale 
temperature 

P 

pressure 
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Tm 

(Tm)^ 

Vs 

P 


pressure  at  the  base  of  a  particular  layer 

the  molecular  scale  temperature 

the  value  of  Tm  at  the  base  of  a  particular  layer 

speed  of  sound  (ft/ sec) 
density  of  the  air  ( slugs/ ft  ) 


The  skeleton  of  the  COESA  is  from  reference  15  and  is  given  in  table  1. 


TABLE  1. 


SKELETON  OF  THE  U.S.  STANDARD  ATMOSPHERE— 1962 


PwflolBf  i«Bp*nitvr«  ttwl  •olMoUr  wlfbit  e»f  tb*  propoM4  0.  8.  Staatori 
«o4  ecaqpaUA  (uA  vter*  •  pKmtrU 

b  •  ctopetMtUl  T  a  blMtle  N  «  mmm  BplMBlak*  wlfbt, 

L  s  inblMi  of  aolMoUr  •e«l«  tooporoturot  »  d1||  /db  (bolev  79  ftop.  to)  ■ 
/tt  (aboT*  79  taop.  to),  7n  s  BoUcuUr  aeaU  t«opoi«t«v*  *  (t^)  M.;  aoi 
Mq  s  valw  of  M* 


*» 

b,  to 

N 

7,  2 

r.  i* 

e.  t/’^ 

0.000 

0.000 

266.15 

-6.5 

28.966 

266.15 

10.1325  /  *• 

1.2250  /  3* 

11.019 

11.000 

216.65 

0.0 

28.966 

216.65 

2.2632  /  2 

3.6392  /  2 

SO.063 

20.000 

216.65 

1.0 

28.966 

216.65 

5.*7*7  /  1 

8.6033  /  1 

32.1C8 

32.000 

228.65 

2.8 

26.966 

226.65 

6.6798  0 

1.3225  /  1 

*7.350 

*7.000 

27D.65 

0.0 

26.966 

270.65 

1.1090  0 

).*275  0 

5*.W9 

92.000 

270.65 

•2.0 

26.966 

270.65 

5.8997  -  1 

7.9939  •  1 

«1.59l 

61.000 

*52.65 

-b.o 

28.966 

252.65 

1.8209  -  1 

2.5106  -  1 

79.99* 

79.000 

lSo.65 

0.0 

28.966 

180.65 

1.0376  -  * 

2.0009  •  2 

90.000 

88.7*3 

160.65 

3.0 

26.966 

180.65 

1.6*37  -  3 

’  1698  -  3 

100.000 

98.*51 

210.65 

5.0 

26.66 

210.02 

3.0070  -  * 

*.9731  -  * 

110.000 

106.129 

260.65 

10.0 

26.56 

297.00 

7.3927  -  9 

9.8277  -  5 

120.000 

117.777 

360.65 

20.0 

28.07 

3*9.*9 

2.5209  -  5 

2. *392  -  5 

190.000 

1*6.5*2 

960.65 

15.0 

26.92 

892.79 

5.0999  -  6 

1.8350  -  6 

ifio.ooo 

156.071 

1,110.65 

10.0 

26.66 

1,062.20 

3.6929  -  6 

1.158*  -  6 

170.000 

165.57* 

1,210.65 

7.0 

26.*0 

1,103.*0 

2.7915  -  6 

8.0330  -  7 

190.000 

18*.*85 

1,350.65 

9.0 

29.85 

1,205.*0 

1.68*5  -  6 

*.3*50  -  7 

230.000 

221.966 

1,550.65 

^.0 

2b.70 

1,322.30 

6.9572  -  7 

1.5631  -  7 

300.000 

266.*76 

1,630.65 

22.66 

1,*3*.10 

1.8628  -  7 

3.5831  -  8 

*00.000 

376.315 

2,160.65 

19.9* 

1,*87.*0 

ii.02T8  -  8 

6.*9*5  -  9 

900.000 

*63.530 

2,*80.65 

17.9* 

1, *99.20 

1.09*9  -  8 

1.5758  -  9 

600.000 

5*6.235 

2,590.65 

1.1 

16.6* 

1,906.10 

3.**75  -  9 

*.69&  -  10 

700.000 

630.536 

2,700.65 

16.17 

1,507.60 

1.1908  -  9 

1.5361  •  10 

•Powr  of  10  bj  toiob  prtetdlng  outoor  wu%  bo  ■nltlpliod. 


Reprinted  with  permission  of  ASTRONAUTICS,  a 
publication  of  the  American  Rocket  Society. 
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This  routine  must  be  called  by  the  CODAP  symbolic  language  as  follows: 


LDA 

ALT 

RTJ 

ATMOS 

TEM 

OCT 

PRES 

OCT 

RHO 

OCT 

VA 

OCT 

ERR 

N,  R. 

where 


ALT 

is 

the  altitude  in  feet 

TEM 

is 

the  temperature 

PRES 

is 

the  pressure 

RHO 

is 

the  density 

VA 

is 

the  speed  of  sound 

ERR 

Is 

the  error  return 

N.R. 

is 

the  normal  return 
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b.  Subroutine  _£  clock  and  subroutine  S  clock. 

These  subroutines  are  incorporated  to  compute  the  time  used  by 
the  computer  in  computing  a  typical  trajectory.  These  subroutines  are 
written  in  the  CODAP  symbolic  language  and  are  callable  by  FORTRAN. 

The  S  clock  subroutine  will  initialize  the  computer  clock  to  zero,  and  the 
E  clock  subroutine  will  stop  the  clock  and  print  out  the  elapsed  time  in  hours, 
minutes,  and  seconds. 

To  use  the  S  clock  subroutine: 

CALL  SCLOCK 

To  use  the  E  clock  subroutine; 

CALL  ECLOCK(N)  where  time  is  printed  on  tape 
number  N. 
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The  subroutine  converts  geocentric  l^ltude  and  radius  to  geodetic 

g 

latitude  and  altitude  by  use  of  the  following  equations. 


where  a^  =  equatorial  radius  of  the  Earth 

f  =  flattening  of  the  Earth 
r  =  geocentric  position  vector 
H  =  geodetic  altitude 
6  =  latitude 

c  =  refers  to  geocentric 
G  =  refers  to  geodetic 

This  subroutine  is  callable  by  FORTRAN. 

To  use: 

CALL  GEODED  (A,  B,  C,  D) 

where  A  Is  the  geocentric  latitude  (radians) 

B  is  the  geocentric  position  vector  (feet) 
C  Is  the  geodetic  latitude  (radians) 

O  is  the  geodetic  altitude  (feet) 
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d.  Impact  subroutine. 

The  Impact  subroutine  is  a  point  mass  three-degree-of-freedom 
trajectory  subroutine.  This  routine  is  incorporated  primarily  to  compute  the 
trajectories  of  the  "separated"  expended  stages. 

The  vector  form  of  the  equation  of  motion  is  the  geocentric  position 
equation  derived  in  section  2  of  this  report  and  is 

d.  R  =  a£.  +  2  *  “it)  - 

dt  M  dt  ^  ^  ^ 

The  atmosphere  subroutine  is  used  along  with  the  oblate  Earth 
described  in  section  2.  The  drag  parameter  (C^A/M)  vs  .  Mach  number 

table  for  the  expended  stages  are  read  in  the  main  program  and  placed  in 
common  for  use  by  the  Impact  subroutine. 

The  Integration  routine  will  use  the  same  option  as  the  main  portion 
of  the  program  (powered  flight)  uses.  When  the  altitude  is  negative,  the 
Impact  subroutine  will  punch  a  card  containing  the  name,  stage  number,  im¬ 
pact  latitude,  and  longitude.  The  routine  will  then  terminate  and  return  to 
the  main  program.  To  use  the  Impact  subroutine,  control  number  8  must 
be  set  equal  to  1  or  to  2. 
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VARIABLES  USED  IN  IMPACT  SUBROUTINE 


AB 

Not  used  In  impact 

AC 

Not  used  in  impact 

ACODE 

integration  code 

ACODES 

Integration  code  table 

AD 

Not  used  in  impact 

AE 

Equatorial  radius  of  the  Earth 

AG 

Not  used  in  impact 

AH 

Not  used  in  impact 

AI 

Not  used  in  impact 

AK 

Not  used  in  impact 

AL 

Block  reserved  for  equivalence 

ALT 

Altitude  of  empty  stage 

AM 

Not  used  in  impact 

AN 

Not  used  in  impact 

ANl 

Variables  used  for  output  parameters 

AN2 

Variables  used  for  output  parameters 

AO 

Not  used  in  impact 

AP 

Not  used  in  impact 

ARC 

Argument  used  for  gravitational  computation 

CBLACK 

Block  reserved  for  integration 

CDM 

Cj^A/M  -  drag  parameter 

CHECK 

Check  used  in  integration 

COSRA 

Cosine  of  the  range  angle 

DEL  PR 

Print  time  increment 

DM 

Drag  over  mass  parameter 

DRA 

Variable  used  in  computing  drag 

DRAG 

Drag  parameter  of  empty  stage 

DUMBl 

Not  used  in  impact 

ERR 

Integration  error 

FMN 

Mach  number  of  empty  stage 

*  Stored  in  common 
f  Equivnlonc*  with  AL 
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GM 

GX 

GY 

GZ 

H 

1 

ICODE 

II 

JP 

J,T 

J6 

K 

KAA 

KAf 

KAZ 

N 

NAME 

NN 

NOE 

NOT 

OUT 

PE 

PI 

PR  TIM  1 
R 

RAD 

RANGE 

R2 

SIT 

SMAX 

SMIN 


Gravitational  constant 
X  component  of  gravitational  attraction 
Y  component  of  gravitational  attraction 
Z  component  of  gravitational  attraction 
Variable  used  for  first  time  Increment 
Utility  index 

Code  for  method  of  integration 
Utility  index 
Page  count 

Number  of  stage  being  computed 

Number  of  lines  per  page 

Utility  index 

Not  .used  in  impact 

Not  used  in  impact 

Not  used  in  impact 

Number  of  values  in  drag  table 

Name  stored  in  common 

Integer  for  drag  selection 

Number  of  equations  for  integration 

Number  of  output  tape 

Dimensio.ied  output  variables 

Earth  flattening  constant 

X  =  3.1415927 

Print  time  for  integration  routine 
Length  of  position  vector 
=  x/180,  =  1/57.2957795 

Great  circle  range  from  launch  point 
Position  vector  squared 
Sine  of  the  latitude 
Maximum  integration  step  size 
Minimum  integration  step  size 


+ 


* 

* 

* 

* 

* 

+ 

* 


* 


+ 

+ 


*  Stored  in  common 
-f  Equivalence  with  AL 
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SS 

T 

TFIMP 

THET 

TRERR 

VA 

VD 

VEL 

W 

WD 

W2 

X 

XALT 

XD 

XDD 

XDIMP 

XDl 

XIMPA 

XK2 

XL  AT 

XMACH 

XMN 

XX 

XXD 

Y 

YD 

YDD 

YDl 

Z 

ZD 

ZDD 

ZDl 

Z2 


Integration  step  size  + 

Time,  independent  variable  + 

Initial  time  for  Impact  computation  * 

Latitude  of  position  vector 

Integration  truncation  error  + 


Velocity  of  sound 
Vector  used  for  matrix  rotation 
Velocity  of  empty  stage 
Rotational  velocity  of  the  Earth 
Vector  used  for  matrix  rotation 
=  W  squared 


Component  of  position  vector  + 

Altitude  from  GEOD  subroutine 

Velocity  vector  component  + 

Second  derivative  of  the  motion  equation  + 

Input  velocity  vector  * 

Velocity  used  in  integration  routine  + 


Input  position  vector 
Oblateness  constant 

Geodetic  latitude  from  GEOD  subroutine 
Mach  number  table  for  drag 
Stage  Mach  number  table 
Vector  for  computing  range 
Velocity  vector 


Component  of  position  vector  + 

Velocity  vector  component  + 

Second  derivative  of  the  motion  equation  + 

Velocity  used  in  integration  routine  + 

Component  of  position  vector  + 

Velocity  vector  component  + 

Second  derivative  of  the  motion  equation  + 

Velocity  used  in  integration  routine  + 

Number  of  lines  per  page  * 


*  Stored  in  common 
-t  Equivalence  with  AL 
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e.  integration. 

A.  IDENTIFICATION 

TITLE:  Numerical  Integration  of  Ordinary  Differential 
Equations  with  Error  Control 

CO-OP  ID;  D2  CODA  NMI  2 

PROGRAMMER:  Roger  Johnson 

DATE;  February  27,  1961 

B.  PURPOSE 

To  integrate  a  set  of  N  simultaneous  first  order 
differential  equations  of  the  form: 

dx. 

—  f.  (t,  Xj.  x^,  i  =  (1,2,  ...N). 

The  user  has  the  option  of  either  a  Runge-Kutta  or 
Adams-Moulton  integration  scheme.  The  integration  step 
size,  ot,  may  be  a  variable  step  size  under  error  control  or 
a  fixed  step  size.  A  print  option  causes  printing  of  the  initial 
conditions  and  various  parameters  before  the  start  of  the  in¬ 
tegration  to  help  define  each  case.  The  user  also  has  an 
option  to  break  into  the  program  at  exact  specified  values  of 
the  independent  variable  t. 

C.  USAGE 

The  index  registers  are  saved.  No  internal  checks  of 
arithmetic  faults  (exponential  overflow,  underflow,  etc.  )  are 
made.  When  the  print  option  is  used,  the  differential  equation 
routine  removes  the  selection  of  interrupt  on  arithmetic  faults 
before  printing  and  does  a  clear  arithmetic  faults  after  print¬ 
ing.  The  differential  equations  routine  should  not  cause  an 
arithmetic  fault  unless  the  x^  's  exceed  the  range  of  the 
floating  point  format. 
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1.  Calling  Sequence  - 

The  calling  sequence  used  to  start  or  restart;  if  parame¬ 
ters  are  changed,  the  integration  of  a  set  of  differential 
equation  is 


Loc 

om. 

B 

M. 

OPN 

s. 

M 

BJbTA 

SLJ 

4 

ADAMS 

0 

P 

DERIV 

BETA+ 1 

0 

0 

TEXIT 

0 

0 

T 

BETA+  2 

0 

0 

DATA 

0 

0 

COMMON 

BETA+  3 

0 

0 

EXIT 

0 

0 

TMIN 

BETA+  4 

ERROR 

RETURN 

2.  Parameters  - 

The  upper  part  of  the  first  instruction  of  the  calling 
sequence  is  a  return  jump  command  SLJ  4  ADAMS  ,  to  the 
first  location  of  differential  equation  routine. 

If  E  is  unequal  to  zero  the  differential  equation  routine 
will  print  the  initial  conditions  and  the  parameters  in  the  DATA 
region  using  the  Generalized  Listable  Output  Routine 
(J5  LMSD  OUTPUT).  If  E.  is  zero  no  printing  takes  place. 

DATA:  starts  a  block  of  3N-t-  8  locations  that  the  user  sets  up 
with  the  parameters  and  variables.  The  location  and  descrip¬ 
tion  of  the  parameters  and  variables  are  as  follows: 

DATA  contains  the  integration  scheme  code  which  is  a  fixed 
point  binary  integer  (binary  point  at  the  far  right)  set  up  by  the 
user. 

a)  If  code  =  0:  The  integration  scheme  is  the  Runge- 
Kutta  mode  with  a  fixed  At. 

b)  If  code  =  1:  The  integration  scheme  is  the  Runge- 
Kutta  predictor- corrector  mode  with  a  variable  At. 
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c)  If  code  =  Z:  The  integration  scheme  Is  the  Adams- 
Moulton  mode  with  a  fixed  At. 

d)  If  code  =  3:  The  integration  scheme  is  the  Adams- 
Moulton  predictor-corrector  mode  with  a  variable  At. 

e)  If  the  code  is  negative  or  any  binary  digit  greater  than 
3,  the  code  is  out  of  range.  If  this  occurs,  the  AC  is 
set  to  the  binary  integer  2  and  a  jump  is  made  to  the 
error  return  in  BETA+  4. 

DATAt  1  contains  the  floating  point  number  A  used  in  the 
variable  step  mode  to  prevent  unnecessary  halting  of  At 
when  becomes  small  in  magnitude.  If  A  is  positive  then 
a  positive  floating  point  number  A^  must -be  determined  for 
each  x^,  (i  =  1,2,...  N),  and  be  set  up  by  the  user  in  the 
following  locations: 

Aj  in  DATA  +  2N  +  8 

A^  in  DATA  +  2N  +  9 
i 

AN  in  DATA  +  3N  +  7 

If  A  is  a  negative  floating  point  number  then  the  absolute 
value  of  A  is  set  in  location  DATA  +  2N  +  8  and  used  for 
all  A.  (i=  1,2,...,  N).  If  A  is  zero  then  location  DATA  + 

2N  +  8  is  set  to  zero  and  A  is  ignored  in  the  halting  and 
doubling  option. 

DATAt  2  contains  a  positive  floating  point  number  E  used  in 

the  truncation  error  test  in  the  predictor-corrector  variable 

At  mode.  If  E  is  set  to  10  approximately  h  significant 

figure  are  asked  for  in  the  truncation  error  test. 

- 1  -8 

(10  <E<10  )  is  the  suggested  range  of  £. 

DATAf  3  contains  a  positive  floating  point  number  called 
MINIMUM  OT.  If  in  the  truncation  error  test,  after  a  step  At, 
one  or  more  of  the  variables  doesn't  meet  the  convergence  test 
and  At  is  less  than  or  equal  to  MINIMUM  DT,  the  differential 
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equation  routine  does  a  return  jump  to  the  users  subroutine 
starting  at  TMIN.  The  location  of  TMIN  Is  given  In  the  calling 
sequence.  The  user  may  stop  and  do  some  checking  or  do  an 
unconditional  jump  to  TMIN  to  return  to  the  differential  equa¬ 
tion  routine.  The  differential  equation  after  the  return  accepts 
the  step  At  Ignoring  the  failure  of  the  convergence  test  and 
continues  the  Integration. 

DATAt  4  contains  the  positive  floating  point  number  MAXIMUM 
DT.  If  In  the  truncation  error  test  all  the  variables  have  con¬ 
verged  so  well  that  doubling  Is  Indicated  but  AT  Is  already 
greater  or  equal  to  MAXIMUM  DT  then  AT  Is  not  doubled  and 
the  next  Integration  step  Is  still  AT. 

DATA-t-  5  contains  a  positive  fixed  point  integer  N, ,  the  number 
of  differential  equations  to  be  solved. 

DATA-t  6  contains  the  floating  point  number  DELTA  T.  In 
the  fixed  AT  mode  the  whole  Integration  Is  done  with  the 
fixed  Integration  step  DELTA  T.  In  the  halflng  and  doubling 
mode  the  Initial  trial  step  Is  AT  but  If  the  convergence  test 
falls,  the  Initial  step  will  be  redone  at  half  AT  and  this  half¬ 
lng  will  continue  until  the  convergence  test  Is  passed.  If  the 
convergence  test  Indicates  doubling  the  next  step  made  will 
be  2 AT.  If  on  entering  the  differential  equation  routine 
DELTA  T  Is  zero  or  negative  the  AC  Is  set  to  the  binary 
integer  1  and  an  unconditional  jump  Is  made  to  the  error 
return  BETA  -1  4. 

DATA-t-  7  Is  the  location  of  Independent  variable  t.  The  user 
must  set  up  an  initial  value  of  t  (a  floating  point  number)  that 
can  be  negative,  zero,  or  positive.  The  differential  equation 
will  advance  t  by  some  At  during  each  Integration  step. 

DATA-t-  8  Is  the  beginning  of  the  N  locations  containing  the 
dependent  variable  x^  through  x^^.  The  user  sets  them  up 
to  their  initial  values  at  the  start  of  a  problem.  The  dif¬ 
ferential  equation  integrates  them  with  step  At  and  replaces 
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them  with  their  new  values. 

is  in  location  DATA  -I-  8 
Is  in  location  DATA  +  9 

Xj,^  is  in  location  DATA  +  N  +  7 

DATA  t  N  8  is  the  beginning  of  the  N  locations  of  the 
derivatives  of  x.  or  L(t,  x^,  x^, ....  x^).  The  user  must  code 
a  subroutine  starting  at  location  DERIV  (given  in  the  calling 
sequence)  to  calculate  the  derivatives  of  using  the  values 
of  t,  Xj.x^,  through  Xj^  in  their  DATA  locations  and  then 
store  the  derivatives  in  their  locetions  in  the  DATA  region. 

As  the  DERIV  subroutine  is  entered  by  a  return  jump  from  the 
differential  equations  routine,  an  unconditional  jump  to  location 
DERIV  gives  the  return  to  the  differential  equations  routine  so 
it  can  continue  the  solution.  The  DATA  locations  of  f.,  or 
the  derivatives  are; 

fj  is  in  location  DATA  +  N  +  8 

f^  is  in  location  DATA  +  N  +  9 

f.,  is  in  location  DATA  +  2N  +  7 
N 

DATA  4-  ZN  t  8  is  the  beginning  of  the  N  locations  of  A^  . 
These  are  the  last  locations  used  in  the  DATA  region.  A 
description  of  the  A^'s  has  already  been  given  under  the 
discussion  of  DATA  -t-  1. 

DERIV  is  the  beginning  of  a  subroutine  to  be  coded  by  the  user. 
The  user's  DERIV  subroutine  function  uses  the  variables  t  and 
x^  in  the  DATA  regions  to  calculate  the  derivatives,  f^  ,  and 
stores  them  In  the  DATA  region.  The  locstion  of  the  variables 


dt 


^  = 
dt 
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in  the  DATA  region  have  been  described  before,  but  to  repeat; 

t  is  in  DATA  +  7 

is  in  DATA  +  8 

is  in  DATA  +  9 

Xj^  is  in  DATA  +  N  +  7 

is  in  DATA  +  N  +  8 

is  in  DATA  +  N  +  9 

£j^  is  in  DATA  +  2N  +  7 

After  calculating  and  storing  the  f^'  s  the  DERIV  routine 
does  an  unconditional  jump  to  DERIV  which  causes  a  return 
to  the  differential  equations  routine  so  it  can  continue  the 
solution.  No  printing  should  be  done  in  this  routine  as  the 
x.'s  in  the  DATA  region  may  be  just  preliminary  estimates. 

EXIT  is  the  beginning  of  a  subroutine  coded  by  the  user  to 
perform  printing.  When  t  and  the  xJs  have  been  integrated 
a  step  At  and  the  x^'s  have  satisfied  the  convergence  condi¬ 
tions,  the  differential  equation  routine  goes  to  the  DERIV  sub¬ 
routine  which  calculates  the  derivatives  of  the  new  x.'s.  The 

1 

differential  equation  routine  then  executes  a  return  jump 
instruction,  SLJ  4  EXIT,  to  the  users  subroutine  for  possible 
printing.  The  user's  subroutine  can  just  bump  a  counter  and 
do  printing  just  every  k  steps  or  can  print  at  each  step.  The 
user  can  take  selected  values  of  t,  x.  and  the  derivatives  of 
x^  from  the  DATA  region,  or  calculate  functions  of  these 
variables  and  print  them  or  save  them  for  future  interpolation. 

If  the  Generalized  Listable  Output  Routine  is  used  the  user 
should  remove  the  selection  of  interupt  on  arithmetic  fault 
before  printing  as  the  output  routine  often  causes  arithmetic 
faults  not  related  to  computational  errors.  After  printing,  a 
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clear  arithmetic  faults  instruction  should  be  given  to  clear 
possible  arithmetic  faults  from  the  output  routine.  At  the  end 
of  the  EXIT  subroutine  the  user  does  an  unconditional  jump  to 
EXIT  to  return  control  to  the  differential  equation  routine. 

TEXIT  and  T  are  set  up  by  the  user  to  get  points  on  the 

solution  of  the  differential  equation  for  specified  values  of  the 

independent  variable  t.  The  way  it  works  is  if  the  independent 

variable  t  is  equal  to  or  greater  than  the  contents  of  location 

T  ,  the  differential  equation  routine  will  do  a  return  jump  to 

TEXIT  with  the  contents  of  location  T  minus  the  independent 

variable  t  in  the  AC.  In  using  this  feature  the  contents  of 

location  T  should  be  always  set  greater  than  the  independent 

variable  t.  Because  if  the  next  integration  step  will  make  t 

larger  than  the  contents  of  location  T  the  differential  equation 

routine  does  a  special  Runge-Kutta  integration  step  with  a 

value  of  At  to  advance  t  and  calculate  x.  and  the  deriva- 

1 

tives  of  X..  At  is  such  that  t  is  equal  to  the  contents  of 
location  T.  After  this  special  At  step  the  old  At  is  re¬ 
stored  and  the  differential  equation  routine  does  a  return 
jump,  SLJ  4  TEXIT,  to  the  user's  TEXIT  subroutine  and  as 
t  equals  the  contents  of  location  T  the  AC  is  zero.  The 
user's  subroutine,  TEXIT,  may  do  printing,  keeping  in  mind 
the  suggestions  given  in  the  EXIT  subroutine.  It  should  then 
advance  the  contents  of  T  to  the  next  break  point.  To  exit, 
an  unconditional  jump  to  TEXIT  returns  to  the  differential 
equation  routine. 

If  the  contents  of  T  are  less  than  the  independent  variable 
t  because  of  improper  updating  or  some  other  oversight,  the 
differential  equation  routine  still  executes  the  instruction, 

SLJ  4  TEXIT,  but  with  the  AC  minus  to  show  that  the  break 
point  has  been  passed  in  t.  The  break  point  features  of  the 
differential  equation  routine  will  be  ignored  if  T  or  the 
contents  of  location  T  are  set  to  zero,  then  no  jump  to 
TEXIT  will  ever  occur.  If  both  T  and  the  contents  of 
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location  T  are  nonzero,  but  TEXIT  is  zero,  the  differential 
equations  routine  will  go  to  its  error  return  with  a  fixed 
integer  -1  in  the  AC  because  no  subroutine  can  start  at  zero. 

If,  when  t  equals  the  contents  of  location  T  the  formulas 
to  compute  the  derivatives  of  x  are  changed  or  it  is  wished  to 
change  some  of  the  parameters  in  the  DATA  region,  the  solu¬ 
tion  of  the  differential  equation  must  be  restarted  by  either 
going  to  a  new  calling  sequence  or  back  to  the  start  of  the  old 
calling  sequence. 

If  t  is  at  the  end  of  a  case,  the  next  case  can  be  set  up  and 
the  differential  equation  routine  restarted  with  a  calling 
sequence.  Only  in  the  user's  EXIT  and  TEXIT  subroutines  can 
the  values  of  the  variables  be  trusted  either  for  printing  or  re¬ 
starting  the  solution  by  going  to  a  calling  sequence. 

BETA  -t  4  is  the  error  return  for  the  differential  equation 
routine.  An  i^nconditional  jump  is  made  to  BETA  -f  4,  with  a 
fixed  point  binary  integer  in  the  AC  which  tells  the  type  of 
error,  if  a  parameter  is  obviously  bad.  To  correct  this,  the 
bad  parameter  should  be  changed  and  the  case  rerun.  When  at 
location  BETA  ■¥  4,  if: 

AC  is  0:  then  N,  the  number  of  equations  to  be  solved, 
(location  DATA  -I-  5)  has  been  set  to  zero,  a  negative 
number,  or  is  greater  than  200. 

AC  is  +  1:  then  At  (location  DATA  +  6)  is  either  zero, 
negative,  or  not  in  normalized  floating  point  form. 

AC  is  2:  then  the  integration  code  in  location  DATA  is 
either  negative  or  greater  than  -t-  3. 

AC  is  -1:  Both  T  and  the  contents  of  T  are  nonzero, 
but  TEXIT  is  zero.  This  means  the  exact  t  break  feature 
is  being  used  but  the  users  subroutine,  TEXIT,  starts  at 
location  zero.  This  is  invalid. 
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COMMON  is  a  block  of  14N  -t-  5  locations  starting  at  COMMON 
that  the  user  must  reserve  for  the  differential  equation  routine. 
When  command  has  been  transferred  to  the  users  EXIT  sub¬ 
routine  the  first  N  location  of  COMMON  contain  the  predictor 
values  of  the  new  .  The  truncation  error  of  the  last  step 

for  the  variable  x.  is  about  1/14  (x  .-x.)  in  magnitude.  The 
i  '  p,  i  i  ® 

locations  of  the  variables  x  .  and  x.  are: 

P.  1  I 

X  ,  is  in  location  COMMON 
P.  1 

X  ,  is  in  location  COMMON  ■¥  1 
P.  2 


X  is  in  location  COMMON  +  N-1 

p,  N 

If  the  user  needs  the  values  of  t  and  the  variables  x.  of  the 

1 

lASt  step  for  interpolation  they  can  be  found  in  the  locations: 

t  in  COMMON  +  N 
x^  in  COMMON  +  N  +  1 

x^  in  COMMON  +  N  +  2 


Xj^  in  COMMON  +  2N 

As  already  described  in  the  writeup,  the  current  x^'s  are  in 
the  DATA  locations: 

x^  in  DATA  +  8 

X2  in  DATA  +  9 


in  DATA  +  N  +  7 
N 

TMIN  is  the  first  location  of  a  subroutine  coded  by  the  user 
that  the  differential  equation  goes  to  if  the  convergence  test 
has  failed  and  At  is  less  than  or  equal  to  MINIMUM  DT. 
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Its  use  is  discussed  later  in  this  writeup. 

3.  SPACE  REQUIRED  -  is  about  700  locations  not  including  the 
610  locations  of  the  Generalized  Listable  Output  routine. 

4.  TEMPORARY  STORAGE  -  none.  The  two  blocks  of  storage 
DATA  (3N+  8  locations)  and  COMMON  (14  N+ 5  locations)  can¬ 
not  be  used  by  the  programer  to  store  numbers. 

7.  ERROR  STOPS  -  the  error  return  is  described  before.  It  is 
in  location  BETA  -t  4  in  the  calling  sequence. 

9.  PRINT  INFORMATION  -  If  the  print  option  is  used  the 

"General  Listable  Output  Routine"  must  be  in  the  computer  and 
its  symbol  location,  OUTPUT,  in  70412B.  At  the  head  of  the 
assembly  is  the  card 

OUTPUT  EQU  70412B 

The  symbolic  location  ADTAPE,  in  the  differential  equation 
routine  determines  the  channel  number  of  the  output  tape,  the 
tape  number,  and  1607  number.  The  print  option  sets  up 
words  for  the  write  BCD  tape  calling  sequence  by  using 
ADTAPE  after  inserting  the  proper  carriage  control  in  the 
lower  address.  Therefore,  to  change  the  channel  number, 
tape  number,  or  1607  number,  just  modify  the  contents  of 
ADTAPE  in  the  differential  equation  routine.  At  present: 

CN  or  channel  number  of  the  output  tape  is  4 

TN  or  tape  number  of  the  output  tape  is  4 

UN  or  1607  number  is  2 


D.  METHOD 

To  compute  the  change  of  with  the  integration  step  At  the 
fourth  order  Runge-Kutta  method  gives  the  formulas: 
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=  ( 

At)  L 

(t. 

*1,X2** 

•V 

=  ( 

At)  f. 

(t 

+  ViAt,  X 

1  ^1,  r  •  ■ 

k_  . 

3, 1 

=  ( 

At)  f. 

(t 

+  l/iAt,  x^ 

+  lA  k^^  ^, . . . 

k^,  i 

=  ( 

At)  f. 

(t 

+  At,  x^ 

^3  ^ 

n  ^  ^3,  n) 

Ax. 

1 

=  I 
6 

<^.i 

+ 

2k  ,  .  + 

2,  i 

2k,  .  +  k.  .) 
3, 1  4,  i 

The  c 

alculation: 

c 

X. 

« 

■  (t  +  At) 

=  x.(t)  +  Ax. 

1  1 

is  done  in  doubb 

i  i  i 


precision  to  help  control  rounding  errors.  All  other  calculations 
are  performed  in  single  precision.  Programers  at  STL  didn't  feel 
the  Gill  or  other  modifications  to  control  rounding  errors  to  be  of 
much  value  so  the  standard  Runge-Kutta  method  was  programed. 

The  Adams-Moulton  method  requires  the  derivatives  of 
three  equal  steps  behind.  To  start  the  Adams-Moulton  integration, 
several  Runge-Kutta  integrations  are  required.  The  Adams- 
Moulton  predictor-corrector  formulas  are: 


'i,  k+  1 


=  X. 

i 


.  +  (55  f.  ,  -  59  f.  ,  ,  +  37f. 

k  24  i,  k  i,  k-1  i, 


1,  7-9  f- 
k-2  i. 


k-3' 


Ax 


i.  k 


Ai 

24 


<9  i?  w  1 

i,  k+  1 


+  19  f.  .  -5  f.  ,  1+  f.  , 

i,  k  1,  k-1  i,  k-2 


Q 

Again  the  calculation  x.  .  .  ,  =  x.  i  .  ,  +  Ax.  ,  is  done  in  double 

X|  X  If  Ki  i  If  K 

precision  to  help  control  rounding  errors. 

Both  Adams-Moulton  and  Runge-Kutta  have  error  terms  of  the 
order  (  At  but  the  Runge-Kutta  step  requires  four  references  to 
the  derivatives  against  just  two  references  in  an  Adams-MOulton 
step;  so  for  most  problems  the  Adams-Moulton  integration  method  is 
to  be  preferred. 

In  integrating  a  differential  equation  numerically,  it  is  re¬ 
placed  with  a  difference  equation  and  we  solve  this  difference  equation. 
If  the  Integration  steps  used  are  small  and  an  integration  scheme  like 
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Adams-Moulton  or  Runge-Kutta  which  have  favorable  stability 
properties  is  used,  the  solution  of  the  difference  equation  is  usually 
close  to  the  solution  of  the  differential  equation.  However,  if  some 
of  the  variables  given  by  the  differential  equation  routine  have 
either  odd  oscillations  or  increase  very  rapidly  with  the  physical 
problem  not  suggesting  such  behavior,  the  solution  is  probably 
unstable  and  greatly  in  error.  In  long  problems  with  many  integra¬ 
tion  steps  a  slight  instability  will  cause  the  solution  to  slowly  drift 
from  the  true  solution  of  the  differential  equation.  An  unstable 
solution  can  often  be  made  stable  by  integrating  with  a  smaller  step 
At  or  by  forcing  a  smaller  step  by  asking  for  more  accuracy  in  the 
predictor-corrector  mode.  The  reason  for  this  is  that  different 
step  sizes  change  the  parameters  relating  to  the  stability  of  the 
difference  equation.  Therefore,  if  two  solutions  with  different  step 
sizes  are  similar  and  close  together,  both  solutions  can  be  accepted 
as  correct.  The  routine  also  makes  it  easy  to  do  checking  by  running 
the  same  case  over  using  both  the  Runge-Kutta  and  Adams-Moulton 
integration  methods.  Usually  the  truncation  error  requirements  in 
solving  a  set  of  differential  equations  dictate  integration  steps 
sufficiently  small  to  insure  stability. 

HALFING  AND  DOUBLING  MODE 


When  doing  an  Adams-Moulton  integration  step  At  to  advance 
the  dependent  variables  from  x^(t)  to  x.(t  +  At)  for  each  variable,  we 
first  calculate  the  predicted  value  x?  (t  -t-  At)  for  each  variable. 

The  users  DERIV  routine  is  used  to  calculate  the  predictor  deriva¬ 
tives,  f.(t  +  At,  x^  (t  +  At),  . . . ,  x^  (t  +  At) ).  Using  the  predicted 
derivatives  the  step  is  finished  by  calculating  (t  +  At)  the 
corrector.  An  estimate  of  the  magnitude  of  the  truncation  error  of 
each  variable  in  the  last  step  is: 


x^  (t  +  At)  -  xP  (t  +  At) 


To  get  an  estimate  of  the  truncation  error  in  the  Runge-Kutta  mode 
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we  first  do  a  Hunge-Kutta  integration  step  of  2At  to  get  the  pre¬ 
dictor,  X?  (t  -fZAt).  Then  restarting  with  the  variables  x^(t)  we  do 
two  Runge-Kutta  integration  steps  each  of  length  At  to  get  the 
corrector,  x.(t  +  2At).  An  estimate  of  the  msgnitude  of  the  trunca¬ 
tion  error  for  each  variable  in  the  step  from  x.(t)  to  Xj^(t  +  2At)  is: 

i —  |x*r  (t  +  2At)  -  X?  (t  +  2At)|  . 

15  I  1  »  1 

Going  from  x^(t)  to  x.it  +  2At)  requires  4  references  to  the 
derivatives  using  the  Adams-Moulton  integration  scheme  and  12 
references  to  the  derivatives  using  the  Runge-Kutta  scheme. 

The  convergence  test  uses  the  parameters  A.  and  E 
described  before  in  the  writeup  of  the  DATA  region.  Using  the  pre¬ 
dictor,  X?  ,  and  corrector,  x^  ,  of  each  variable  from  the  last 
integration  step  (for  either  the  Runge-Kutta  or  Adams-Moulton 
method)  if  the  inequality: 


c 

X  . 
_ L 

- 

xP 

_ 

max 

A., 

xP 

*i 

• 

1 

is  satisfied  for  all  i  (  i  =  1,  2,  . . . ,  N)  then  we  say  the  convergence 
test  has  been  passed  and  the  results  of  the  last  integration  step  will 
be  accepted  by  the  differential  equation  routine.  If  the  inequality 
doesn't  hold  for  any  one  of  the  i  ,  then  the  convergence  test  has 
failed  and  the  results  of  the  last  integration  step  are  not  accepted. 

If  E  in  the  above  inequality  is  replaced  by  E/M,  "where  M  is  in 
symbolic  location  ADC  -t-  2  and  set  to  the  value  20  in  floating  point 
format",  and  the  above  inequality  is  still  satisfied  for  every  i,  we 
say  the  convergence  test  to  double  has  passed.  Again,  if  when  E  is 
replaced  by  E/M  in  the  above  inequality  and  the  inequality  doesn't 
hold  for  any  one  of  the  i  we  say  the  convergence  test  to  double  has 
failed. 

If  the  convergence  test  is  satisfied  but  the  convergence  test  to 
double  has  failed,  then  the  integration  step  At  is  left  unchanged  and 
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the  differential  equation  routine  goes  to  users  subroutine  EXIT  for 
possible  printing.  Upon  return  to  the  differential  equation  routine 
it  integrates  the  variables  again  by  the  same  step  At. 

If  the  convergence  test  to  double  is  satisfied  the  differential 
routine  goes  to  the  users  subroutine  EXIT  for  possible  printing  of 
the  accepted  but  on  return  a  set  up  may  be  made  to  make  the 
next  integration  step  2At.  Even  though  the  convergence  test  to 
double  is  satisfied,  sometimes  At  is  not  doubled.  For  instance, 
if  At  is  already  greater  than  or  equal  to  the  number,  MAXIMUM  DT 
in  DATA  -t  4,  then  At  is  not  doubled.  If  At  was  halfed  within  the 
last  four  steps  in  the  Runge-Kutta  mode,  controlled  by  the  fixed  point 
binary  number  4  in  symbolic  location  ADC,  no  doubling  of  At  is 
done  or  if  At  was  halfed  within  the  last  six  steps  in  the  Adams- 
Moulton  mode,  controlled  by  the  fixed  point  binary  number  6  in  sym¬ 
bolic  location  ADC  +  1,  no  doubling  is  done.  The  above  delays  are 
inserted  to  save  machine  time  that  might  be  wasted  in  halfing  and 
doubling  oscii’ »i.tions. 

In  the  Runge-Kutta  mode  doubling  can  take  place  every  inte¬ 
gration  step  but  in  the  Adams-Moulton  mode  sufficient  delays  may 
not  exist  to  do  the  next  integration  step  at  ZAt.  To  integrate  with 
the  next  step  ZAt;  first  x.(t  +  At)  becomes  x^(t*),  where 
t*  =  t  +  At  the  advanced  independent  variable,  and  x.  (t-5At) 
becomes  }t^(t^-3(ZAt) ).  Therefore,  after  any  halfing  or  doubling 
operation  in  the  Adams-Moulton  mode  At  cannot  be  doubled  until 
3  integration  steps  of  the  same  length  are  made  to  get  the  delays 
needed  for  the  next  integration  step  ZAt. 

If  the  convergence  test  has  failed  At  is  usually  halfed.  The 
only  case  where  At  is  not  halfed,  is  if  At  is  less  than  or  equal  to 
the  floating  point  number,  MINIMUM  DT  set  up  by  the  user  in  location 
DATA  -t-  3.  As  described  before  in  this  writeup  then  At  is  supposed 
to  half  but  it  is  already  less  than  or  equal  to  MINIMUM  DT,.  If  a  return 
jump  is  made  to  the  users  subroutine  TMIN  and  if  the  user  returns  to 
the  differentisl  equation  routine  by  an  unconditional  jump  to  TMIN, 
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then  x^(t  -t  At)  is  accepted  and  the  next  integration  step  will  also  be 
At. 

If  the  convergence  test  has  failed  and  At  is  greater  than 
MINIMUM  DT  the  x.(t  +  At)  in  the  Adams-Moulton  mode  and 
x^  ( t  +  2  At)  in  the  Runge-Kutta  mode  are  not  accepted  and  the 
differential  equation  routine  goes  back  to  variables  x^(t)  and  the 
next  integration  step  is  .  5  At.  In  the  Runge-Kutta  mode  the  differ¬ 
ential  equation  routine  restores  the  old  Xj^(t  )  and  its  derivatives 
it  had  saved  in  COMMON  storage  and  takes  the  saved  result  of  the 
first  RK  step  x^(t  +  At)  and  makes  it  the  new  predictor.  It  then 
enters  the  RK  predictor- corrector  mode  at  the  point  it  computes 
the  two  Runge-Kutta  steps,  now  each  of  length  .5  At  to  get  the 
corrector. 

In  the  Adams-Moulton  mode  after  At  is  halfed  the  derivatives 
of  x^(t  -  .5  At)  and  x^lt  -  1,5  At)  are  needed  for  the  next  integration 
step  of  length  .5At.  Using  Xj(t),  x.(t  -At),  x^<t-2At),  x^(t-3At),  and 
x^(t-4At)  in  the  five  point  interpolation  formula  of  Lagrange; 

x.(t  -.5At)  and  x.(t-1.5At)  are  computed.  The  five  point  formula 

^  5 

was  used  as  its  error  term  is  of  the  same  order  (At)  as  the  error 

terms  in  the  Runge-Kutta  and  Adams-Moulton  integration  methods 

with  integration  steps  of  length  At.  After  interpolation,  the  routine 

restores  x^(t)  and  returns  to  the  Adams-Moulton  mode  in  the 

routine  to  the  part  where  the  Adams-Moulton  step  will  be  done  over 

with  the  new  step  .5At.  As  referred  to  before,  no  doubling  of  At 

will  take  place  after  a  halting  in  Adams-Moulton  mode  until  sufficient 

delays  exist  for  doubling  At.  Tags  were  set  to  delay  4  steps  in  the 

Runge-Kutta  mode  and  6  steps  in  the  Adams-Moulton  mode  to  save 

machine  time. 

The  five  point  Lagrange  Interpolation  Formulas  came  from 
page  118  of  "Introduction  to  Numerical  Analysis"  by  F.  Hildebrand. 

In  the  same  book  can  be  found  the  formulas  of  Runge-Kutta  on  page 
237  and  Adams-Moulton  on  page  200. 
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SETTING  PARAMETERS 

In  solving  a  set  of  differential  equations^we  are  interested  in 
keeping  the  error  at  all  points  of  the  solution  within  certain  specified 
limits.  However,  by  setting  the  parameters  A.  and  E  ,  the 
truncation  error  is  only  controlled  for  each  single  step  in  the  pre¬ 
dictor-corrector  mode.  The  total  error  is  a  combination  of  the 
truncation  and  rounding  errors  of  many  steps.  The  truncation 
error  per  step  for  the  variable  x.  is  less  then  ^  Ex.  if  the 

absolute  value  of  x^  is  greater  than  the  value  of  or  less  than 
^  E  A^  if  Aj^  is  larger  of  the  two.  As  most  problems  are  only  a 

few  thousand  integration  steps  and  are  well  behaved,  the  truncation 
error  per  step  is  suggested  to  be  set  to  one  fiftieth  of  the  total  error 
allowed. 


Two  examples  for  setting  A  and  E  are: 


1)  Integrate  x  =cos  t  with  truncation  error  per  step  less 
-5 

than  10  .  As  the  truncation  error  requirement  is  independent  of 

the  magnitude  of  the  variable  x  ,  set  A  to  1  which  is  equal  to  or 
greater  than  the  absolute  value  of  x  at  any  time. 


As  MAX  (A, 


c 

X. 

I 


xP 

I 


)  =  1.  set  E  =  14  ( 10“^). 


Z)  Integrate  x  =  e*  ,  where  x(o)  =1.,  with  truncation  error 
per  step  less  than  x(10  ^).  As  the  truncation  error  requirement  per 
step  is  some  proportion  of  the  variable  x  set  A  to  zero.  As 
MAX  (A,  I I  t  I  1  ^  ~  ^  "  14(10”^).  If  when  Integra- 

ting  the  last  case  to  t,  seconds  we  set  A  =  e  as 

I  ^  I  t  ^  1 1 


MAX  (A,  X 


A  =  e  we  would  have  large  ^t 


steps  at  the  beginning  resulting  in  great  relative  errors  that  would 
be  carried  through  the  whole  solution. 


To  avoid  unnecessary  halting  of  At;  each  value  of  A^  should 
be  set  slightly  less  than  the  average  magnitude  of  its  corresponding 
variable  x^,  but,  if  the  solution  tends  to  be  inaccurate  or  unstable 
when  some  of  the  variebles  are  small  the  value  of  A^  corresponding 
to  these  variables  must  be  decreased. 
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(  ) 


When  trying  to  see  the  effects  of  the  range  of  a  parameter  to 

check  for  the  total  error  of  an  integration;  change  the  parameter  at 

least  by  a  factor  of  10  to  make  sure  it  reruns  with  a  different  step 

size  and  change  E  or  some  of  the  only  to  isolate  the  effect. 

When  running  a  problem  the  user  is  not  too  familiar  with;  a  typical 

case,  or  the  extreme  cases  plus  a  few  middle  cases,  should  be  run 

with  different  values  of  the  parameters  and  E.  Then  after 

comparing  the  solutions  of  the  same  case  which  has  been  run  with 

different  parameters  A.  and  E,  run  the  rest  of  the  cases  with  the 

values  of  the  parameters  A.  and  £  that  gave  satisfactory  total 

error  bounds.  Sometimes  when  E  is  decreased  the  greater  number 

of  integration  steps  can  increase  the  total  rounding  error  until  it  is 

so  much  greater  than  the  total  truncation  error  that  the  solution  with 

the  smaller  E  has  greater  total  error  than  the  solution  with  the 

larger  E;  however,  as  the  variables  x^  are  accumulated  in  double 

precision  and  since  the  Control  Data  computer  has  a  36  bit  mantissa 

-  8 

then  only  when  E  becomes  less  than  10  are  the  total  errors 
likely  to  increase  when  E  is  further  decreased. 

When  solving  some  problems  a  part  of  the  solution  may  have 
very  small  truncation  errors;  then  At  can  become  so  large  that  the 
solution  becomes  unstable.  The  parameter  MAXIMUM  DT,  in 
DATA  +  4,  is  used  to  prevent  At  from  becoming  too  large.  As  an 
example,  when  a  trajectory  was  being  calculated  in  a  region  of  very 
thin  air,  integration  steps  of  15  seconds  passed  the  truncation  error 
test  but  the  solution  soon  became  unstable  and  did  false  large 
oscillations.  By  setting  a  limit  on  the  integration  step  to  a  few 
seconds  the  rerun  solution  was  stable. 

BREAK  POINTS 

As  described  before  in  the  writeup  the  contents  of  T  and 
TEXIT  can  be  set  up  by  the  user  to  get  points  of  the  solution  of  the 
differential  equation  for  specified  values  of  the  independent 
variable  t. 
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One  use  of  this  feature  is  to  do  printing  every  h  seconds 
of  t.  Initially,  the  contents  of  T  are  set  to  h.  Then  the  k-th 
time  the  differential  equation  routine  does  a  return  jump  to  the 
user's  TEXIT  subroutine  the  user  prints  the  solution,  where  t  =  kh, 
then  bumps  location  T  by  h  and  does  a  return  jump  to  the  differ¬ 
ential  equation  routine  so  it  can  calculate  to  the  next  t  break, 
t  =  (k+l)h. 


Another  use  of  this  feature  is  when  starting  a  solution,  if 
some  of  the  variables  are  Initially  equal  to  zero,  a  different  set  of 
parameters  is  required  for  a  small  value  of  t.  At  the  start  of  the 
problem  the  contents  of  T  are  set  to  a  value  t^  ,  slightly  greater 
than  the  initial  value  of  t  and  the  starting  parameters  are  set  up. 
When  t  equals  t^  the  users  TEXIT  subroutine  is  used  to  modify 
the  parameters. 


When  just  the  contents  of  T  or  any  of  the  convergence 
parameters  like  E,  A.  ,  MINIMUM  DT,  or  MAXIMUM  DT  are 
changed  in  the  middle  of  a  case  by  a  user's  subroutine  the  users 
normal  return,  which  is  an  unconditional  jump  to  the  beginning  of 
his  subroutine  may  be  used  in  returning  to  the  differential  equation 
routine.  However,  almost  any  other  kind  of  change  in  a  users  sub¬ 
routine  requires  a  jump  back  to  the  beginning  of  the  calling  sequence 
or  a  new  calling  sequence.  Such  changes  by  the  user  requiring  a 
return  to  a  calling  sequence  are;  if  the  calling  sequence  is  to  be 
modified,  either  N,  code,  or  At  is  changed,  changes  of  the 
variable  t  or  x^.  This  is  necessary  because  parameter  setups 
from  the  calling  sequence  are  done  only  at  the  beginning  of  the 
routine  and  the  logic  of  the  differential  equation  routine  depends  on 
the  code  to  tell  where  it  has  stored  the  variables  and  delays  of  the 
derivatives  of  x^  in  COMMON.  The  delays  are  also  assumed  to  be 
integral  multiples  of  the  present  At  in  DATA  -t-  6. 

To  return  to  examples  in  the  use  of  the  time  interrupt  feature. 
One  way  to  start  solutions  where  some  of  the  initial  values  of  the 
variables  x^  are  zero  is  to  run  in  the  fixed  At  mode  with  a  small 
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At  until  the  solution  is  started.  Then,  after  the  solution  is 
assumed  started  at  time  t^  change  to  the  predictor-corrector  mode 
by  modifying  the  code,  aid  restart  the  solution  at  ^  entering 

a  new  calling  sequence.  At  the  same  time  TEXIT  could  be  changed 
to  a  new  TEXIT  subroutine  that  does  printing.  To  accurately  simu¬ 
late  a  missile  in  the  velocity  range  of  zero  to  a  few  hundred  feet  per 
second  when  it  is  under  the  influence  of  torques  that  change  its 
direction  (like  crosswinds  or  controllers)  requires  small  integration 
steps  of  the  order  .  01  second  or  large  errors  in  the  orientation  of 
the  missile  occur  from  the  integration  process. 


The  last  use  of  time  interrupt  discussed  is  when  there  are 
discontinuities  of  the  variables  or  their  derivatives  at  specified 
times.  An  example  is  when  a  rocket  engine  burns  out  at  tp.  By 
setting  tp  in  location  T  then  the  user  in  his  TEXIT  routine  must 
modify  his  DERIV  routine  to  omit  the  thrust  term  and  the  change 
the  drag  calculations.  As  such  changes  cause  discontinuities  in  the 
derivatives  of  the  differential  equation  routine  should  be  re¬ 
started  by  going  to  a  new  calling  sequence. 

Another  case  similar  to  the  above  is  in  the  staging  of  missiles, 
there  should  be  a  break  when  a  missile  separates  into  two  sections, 
a  break  at  the  end  of  free  flight  of  the  second  stage,  and  at  the  time 
of  burnout  of  the  rocket  engine;  whose  firing  terminated  the  free 
flight  of  the  second  stage.  When  controllers  are  turned  on  and  off 
they  can  cause  discontinuities  in  the  derivatives  of  the  variables 
x^  or  even  in  x^  if  the  controllers  are  impulse  functions,  so  the 
differential  equation  routine  should  be  restarted  by  entering  a  calling 
sequence  at  such  times.  When  interpolating  through  nonsmooth 
functions  with  discontinuous  derivatives  to  maintain  accuracy  the 
differential  equation  routine  half  s  At  many  times  using  a  great  deal 
of  machine  time;  it  may  also  caise  an  exit  to  the  users  TMIN  sub¬ 
routine  because  At  is  less  than  or  eqial  to  MINIMUM  DT  and  the 
error  requirements  are  not  satisfied.  In  runs  with  controllers  even 
though  the  errors  in  interpolating  variables  through  nonsmooth 
functions  may  be  small  in  terms  of  the  magnitudes  of  the  variables 
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X.,  as  the  error  term  driving  the  controller  is  a  function  of  the 
difference  of  some  ideal  x.  and  actual  x^,  this  error  may  jump 
several  hundred  percent  making  the  rest  of  the  controller  or  guidance 
simulation  useless. 

USERS  INTERPOLATIONS 

The  user  may  require  other  breaks  than  time  breaks.  He  may 
wish  to  print  the  output  every  1000  feet  of  altitude  or  turn  on  or  off  a 
controller  when  the  error,  which  is  a  function  of  the  x.'s,  becomes 

i 

equal  to  some  critical  value.  To  do  this  the  user  could  save  several 
values  of  the  functions  at  different  times  and  for  some  interpolation 
formulas  also  the  derivatives  of  the  functions  and  do  an  inverse 
interpolation  to  find  the  time  when  the  function  of  the  error  reached 
the  critical  point.  The  interpolation  formula  used  should  have  an 

5 

error  term  of  the  order  (At)  to  have  about  the  same  accuracy  as 
the  differential  equation  routine. 

After  the  time  of  the  break  t  (break),  is  known  by  inverse 
interpolation  the  user  can  interpolate  for  all  the  other  variables,  if 
he  has  saved  enough  delays,  then  change  the  method  of  calculating 
the  derivatives  to  include  the  new  status  of  the  controller  and  re¬ 
start  the  differential  equation  routine  by  going  to  a  calling  sequence. 
If  the  user  hasn't  saved  enough  delays  of  the  variables  to  do  suf¬ 
ficiently  accurate  interpolation,  he  could  let  the  differential  equation 
routine  do  the  interpolation.  The  user  is  in  his  EXIT  routine,  the 
time  of  break  is  between  the  current  t  in  DATA  +  7  and  the  t  that 
was  in  DATA  +  7  at  the  time  of  the  last  entry  into  the  EXIT  routine 
as  the  break  occurred  between  steps.  The  user  first  moves  the  old 
t  and  the  old  variable  x.  from  COMMON  in  this  way: 

old  t  from  COMMON  +  N  to  t  in  DATA  +  7 

old  Xj  from  COMMON  +  N  +  1  to  Xj^  in  DATA  +  8 

old  X  from  COMMON  +  2N  to  x  in  DATA  +  N  +  7 
n  n 
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Then  set  the  new  AT  =  t( BREAK)  >  t(old)  in  DATA  +  6  and 
change  the  code  to  zero,  by  setting  zero  in  location  DATA,  and  go 
to  a  calling  sequence  to  restart  the  differential  equations  routine. 

The  first  return  to  users  EXIT  will  be  with  t  =  t  (old)  but  in  the 
second  return  to  EXIT  the  variables  and  x.  will  be  advanced  to  the 

i 

break  point.  The  user  then  makes  the  necessary  changes  in  his 
DERIV  subroutine  that  happen  at  the  break  point,  restores  the  old 
code  and  At  and  goes  to  a  calling  sequence  to  restart  the  differential 
equation. 

The  user  shouldn't  turn  on  the  controller  or  ignite  the  rocket 
motor  in  the  simulations  until  he  has  all  his  variables  interpolated 
at  the  break  time  or  he  will  have  serious  errors  in  his  interpolations 
because  of  the  nonsmooth  functions  introduced.  If  the  user  hasn't 
saved  delays  of  the  function  that  determines  the  break  point  he  could 

calculate  the  old  value  of  the  function  from  t  and  the  variables  x. 

1 

in  COMMON  +  N  through  COMMON  +  2N  and  do  a  linear  interpolation 
for  the  break  point.  Then  advance  t  and  the  variables  to  the 
estimated  break  from  linear  interpolation  by  a  Runge-Kutta  step  as 
before.  Then  the  estimated  t  break  from  linear  interpolation  and 
the  new  variable  x.  with  a  new  interpolation  can  be  used  to  get  a 
better  estimate  of  the  break  time  and  this  process  repeated  until  the 
break  point  is  calculated  to  sufficient  accuracy.  The  user  should 
note  the  process  will  iterate  closer  and  closer  to  the  break  point  but 
may  not  go  past  it  so  the  exit  from  above  process  should  be  done 
when  an  estimate  of  the  t  (break)  is  close  to  the  actual  t  (break). 

TMIN  SUBROUTINE 

If  in  the  predictor- corrector  mode  one  of  the  variables  has 
failed  the  convergence  test  and  AT  is  less  than  or  equal  to 
MINIMUM  DT  the  differential  equation  routine  does  a  return  jump 
to  the  users  subroutine  TMIN. 

If  the  user  has  a  problem  he  knows  can  be  integrated  to 
sufficient  accuracy  using  a  step,  AT,  he  can  set  MINIMUM  DT  to 
AT  and  in  the  TMIN  subroutine  do  an  unconditional  jump  to  TMIN 
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to  cause  the  return  to  the  differential  equation.  The  differential 
equation  routine  accepts  the  last  step  integrated  with  the  last  AT 
and  also  does  the  next  step  with  the  same  AT.  Therefore,  machine 
time  is  not  wasted  with  a  AT  much  smaller  than  AT  and  if  in  other 
regions  of  the  problem,  where  truncation  error  is  small  AT  may  be 
much  larger  than  AT.  The  regions  of  small  AT  could  be  caused  by 
slight  discontinuities  in  curve  fits  of  the  empirical  data  that  determine 
the  coefficients  of  the  differential  equation  or  the  values  of  are 
too  small. 

If  the  user  has  a  problem  that  isn't  too  familiar  MINIMUM  DT 
should  be  set  very  small.  Then  the  TMIN  subroutine  should  stop  the 
problem  so  it  can  be  examined.  The  user  may  find  derivatives 
changing  rapidly  or  which  are  very  large  because  they  are  not  calcu¬ 
lated  correctly  or  an  unnoticed  singular  ‘  is  making  a  derivative 
infinite.  Also  the  user  may  be  changing  t  tables  because  of 

program  errors,  or  if  the  differential  equation  has  impulse  functions 
the  variables  were  changed  but  the  differential  equation  routine 
wasn't  restarted.  If  the  value  of  MINIMUM  DT  was  only  a  factor  of 
10  smaller  than  the  AT  that  runs  most  of  the  problem  accurately,  an 
exit  to  TMIN  could  mean,  discontinuous  derivatives  being  integrated 
through  a  point  where  the  differential  equation  routine  should  have 
been  restarted.  Another  possibility  with  MINIMUM  DT  a  factor  of  10 
smaller  than  the  standard  DT,  if  in  TMIN  when  some  of  the  variables 
are  small;  the  values  of  some  of  the  A^  could  be  increased  and  the 
problem  reran  with  both  a  smaller  MINIMUM  DT  and  the  old  AJ  s 
and  the  old  MINIMUM  DT  and  the  larger  A^'s  and  the  solution  com¬ 
pared  to  see  if  the  shorter  machine  time  in  running  with  larger 
A^'s  gives  sufficiently  accurate  solutions. 

RUNNING  A  PROBLEM  IN  REVERSE 

If  a  user  wishes  to  run  a  problem  backwards  in  time,  as  At 
cannot  be  made  negative  in  this  routine,  the  method  described  below 
must  be  used. 
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Given  the  standard  set  of  N  simultaneous  differential 
equations  to  be  integrated  we  have: 

^  =  f.  Xj(t),  x^lt),  .  . .  Xj^(t)j 

Let  the  initial  value  of  t  be  t  and  its  final  value  be  t.  .  To 

o  f 

run  this  same  set  of  differential  equations  backwards  we  make  the 

transformation  T  =  t^  -  t,  to  the  new  independent  variable  T  . 

Let  the  initial  value  of  T  be  0,  then  t  =  t^  and  the  initial  values 

of  the  variables  x.  is  x.itj).  Then  the  final  value  of  T  is 

t,-t  and  this  corresponds  to  t  =  t  and  the  final  value  of  the 
f  o  ^  o 

variables  x.  is  x.(t  ),  The  derivatives  of  the  variables  x.  with 
i  •  i  o  i 

respect  to  T  are: 

^  [^f  ■ 

as  t  =  tj  -  T  and  dT  =  -  dt. 


E.  TEST  CASE 

To  help  clarify  the  writeup,  a  simple  test  case  is  inserted. 
I  am  integrating  the  three  equations: 


dt 


-4x 


1 


2x 


3 
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Let  Xj^(O)  =  0,  x^tO)  =  Z,  and  x^(0)  =  1.  The  solutions  of  the 
above  differential  equations  are: 

x^  =  sin  2t 

x^  =  2  cos  2t 

2t 

x^  =  e 

Using  the  T  break  feature  printing  is  done  every  second. 

At  t  equal  one  second. 

Xj^  =  sin  2  =  .909297427  the  routine  gives  .9092974248 

x^  =  2  cos2  =  -.832293674  the  routine  gives  .  8322936859 

Xj  =  e^*  =  7.38905610  the  routine  gives  7.389056142 

The  last  printout  is  t  =  5  seconds. 

Xj  =  sin  10  =  -.544021111  the  routine  gives  -.544021143 

x^  =  2  cos  10  =  -1.678143058  the  routine  gives  1.678143027 

Xj  =  e^®  =  22026.4658  the  routine  gives  22026.46649 
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TEST  CASE  CODING 


ORG 

10000 

OUTPUT 

EQU 

70412B 

ADAMS 

EQU 

3720B 

TEST 

ONE 

STA 

T 

FIRST  T  BREAK  IS  1. 

STA 

DATA+  10 

X3  (0)  IS  1. 

LDA 

TWO 

STA 

DATA+  9 

X2(0)  IS  2. 

ENA 

0 

STA 

DATA+  7 

+  0  IS  ZERO 

STA 

DATA+8 

X1(0)  IS  1. 

BETA 

SLJ 

4  ADAMS 

CALLING  SEQUENCE 

ZRO 

1  DERIV 

B  NOT  ZERO  SO  PRINT 

ZRO 

TEXIT 

INITIAL  CONDITIONS 

ZRO 

T 

+ 

ZRO 

DATA 

ZRO 

COMMON 

+ 

ZRO 

EXIT 

ZRO 

TMIN 

A2 

SLS 

A2 

ERROR  RETURN 

ONE 

DEC 

1. 

1. 

TWO 

DEC 

2, 

2. 

SIX 

DEC 

6. 

7. 

MFOUR 

DEC 

-4. 

-4. 

DATA 

DEC 

3 

CODE  IS  3 

DEC 

1, 

A  IN  DATA+  1  IS  PLUS 

DEC 

1.D.8 

E  IN  DATAf  2  IS  .  000001 

DEC 

3.D.4 

MIN  DT  IN  DATA+  3  IS  .  0003 

DEC 

1. 

MAX  DT  IN  DATA+  4  IS  1. 
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3 

N  IS  3  IN  DATA+  5 

DEC 

5.D-3 

INITIAL  DT  IS  .  005 

DEC 

0.  0»  0.  0.  0.  0.  0 

+  XI.  X2.  X3.  and  DERIV  XI.  X2.X3 

DEC 

.1.  .1. 1. 

Al.  A2,  A3 

COMMON 

BSS 

47 

14  W+  5  IS  47 

T 

BSS 

1 

T  BREAK 

TMIN 

SLS 

TMIN 

STOP  IF  DT  TOO  SMALL 

COUNT 

BSS 

1 

DERIV 

SLJ 

0 

LDA 

DATA+  9 

STA 

DATA+ 1 1 

DERIV  XI  IS  X2 

LDA 

MFOUR 

FMU 

DATA+  8 

STA 

DATA+  12 

DERIV  X2  IS  -4X1 

LDA 

DATA+  10 

FAD 

DATA+  10 

STA 

DATA+  13 

2X3  IS  DERIV  X3 

_SIiJ_ 

_ DERIV _ 

EXIT 

SLJ 

0 

EXIT  EACH  DT  STEP 

SLJ 

EXIT 

TEXIT 

SLJ 

0 

SLJ 

4  OUTPUT 

PRINT  TITLE 

A3 

SLS 

A3 

+ 

03 

BCD 

SUBSCRIPT 

■if 

2021 

+ 

03 

BCD-^-2 

Bf 

1029 

X _ 

+ 

03 

BCD+3 

wm 

1047 

DERIV  X 

+ 

03 

BCD-t-4 

n 

1070 

T 

+ 

4  4 

00 

2  10 

PRINT  TITLE 
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ENA 

1 

STA 

COUNT 

ENI 

1 

0 

LOOP 

SLJ 

4 

OUTPUT 

A4 

SLS 

A4 

+ 

04 

COUNT 

COUNT  1  TO  3 

00 

10 

+ 

06 

1 

DATA+  8 

00 

10030 

XI 

+ 

06 

1 

DATA+  1 1 

00 

10050 

DERIV  XI 

+ 

06 

DATA+  7 

00 

10070 

TIME 

+ 

01 

4 

4 

00 

2 

16 

PRINT  X.  DERIV 

X.  AND  T 

RAO 

COUNT 

EVERY  SEC 

+ 

ISK 

1 

2 

BUMP  1 

SLJ _ LOOP _ _ _ 

LDA  T  END  PRINT 


FAD 

ONE 

STA 

T 

FSB 

SIX 

AJP 

M  TEXIT 

STOP 

SLS 

STOP 

IFDATA+7  OR  T  IS  7  STOP 

BCD 

BCD 

2SUBSCRIPT 

_ IX _ 

BCD  1  DERIV  X 


BCD 

1  T 

END 

TEST 
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«  o  o  o 

Q  o  o  o 

-«  o  o  o 

^  ,  o  o  o 

n  +  .000 

■<  <000 

20  000 

^  000 

.  000 


+  +  + 

000 
000 
000 
000 
.000 
H  o  o  o 
000 
000 
000 


+  +  + 

000 
000 
000 
000 
.  000 
H  o  o  o 
000 
000 
000 
rvi  .->]  (M 


000 

000 

000 

000 

000 

Fi  O  O  O 

000 

000 

000 

fO  fO  M 


000 
000 
000 
000 
000 
Fi  O  O  O 

000 

000 

000 

^  ^  ^ 


2  o  O 

3”>2 


0  0^0 

0  nj 

0  CO 

00^ 

+  + 

+  +  + 

+  +  + 

+  +  + 

+  +  + 

0  0 

00  9^  (M 

00  CO  •-! 

rg  0  CM 

0  ID  0 

0  0 

^  in  ^ 

m  fvj  r* 

00  0^  ^ 

nO  fO  >0 

0  0 

ra  00  ^ 

0  rg  0 

in  0 

O'  (M  0 

0  0 

^  >0  >0 

in  h-  in 

^  0  00 

0>  -H  00 

0  0 

X  rg  00  fh 

^  in  ^  00 

..  -H  0  lf> 

0  0 

O'  0^  0 

^  0  rg  00 

fO  eg 

«  ID  0  O' 

0  0 

nj  M  O' 

00  N  9^ 

^  0  ^ 

^00 

0  0 

O'  nj  00 

>9  0  in 

O'  rg  lO 

in  00 

0  0 

0  m  CO 

in  fO  ^ 

e-  0^  0 

00  O'  {^ 

fM  ^ 

00 

n*  in 

CM  ^ 

O'  m  pj 

•  • 

•  •  • 

*  •  • 

•  •  • 

•  •  • 

Pi 

y  -H  rsi « 
M 

Q  (Q 

O 

U 


U  ^  N  rO 

M 

(Q 

3 

V) 


y  -  tv. 

w 

A 

D 

M 


y  ^  M  fo 
OT 

(Q 

3 

(O 


U  I-H  *M  fO 
M 
PQ 
3 

CO 


5440211430  +0  -.1678143115  +  1  .5000000000 

1678143027  +1  .2176084290  +  1  .5000000000 

2202646649  +5  .4405293659  +5  .5000000000 
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f.  i£L££E£F. 

A.  IDENTIFICATION 


TITLE: 

CATEGORY: 

PROGRAMER: 

MODIFIED: 


Divided  Difference  Interpolation  or 
Extrapolation  Routine 

Mathematical  Subroutine 

Sanford  Elkin,  William  Silverman,  and 
Ed  Fleming 

June  1961 


B.  PURPOSE: 


Given  a  table  of  M  values  in  floating  point  which  define  a 
function,  and  given  a  floating  point  argument  X  ,  this  routine 
approximated  the  functional  value  Y  of  that  argument  by  a  poly-  , 
nomial  Interpolation  or  extrapolation. 


C.  USAGE: 


1.  Calling  sequence: 

CALL  INTER PF  (X,  M,  N,  a,  b) 

Interpolate,  find  Y  for  X,  X  Fit  argument 
where  various  X  vs.  Y  given 

in  table.  M  Size  of  the  table 

N  Degree  of  interpolation 
a  Table  of  X  values 
b  Table  of  Y  values 

Where  M  is  in  fixed  integer  form,  N  is  the  order  of  interpolation 
or  extrapolation  desired  —  in  fixed  Integer  form. 

The  X's  and  Y's  must  be  in  floating  format  and  the  table  of  X 
values  must  be  stored  in  decreasing  order. 


105 


TDR-63-11 


D.  MATHEMATICAL  METHOD:  (reference  19) 


Y  =  L,  1,  . . . ,  n 


(x)  _  1 


X  -  x_ 
n  0 


Vl . 


I,  -  (x)  X  -  X 


Y  =  I.  .  (x)  Is  equivalent  to  Lagrange's  Formula 

Y  =  V  (x)  =  (x-  xi)(x-  x^)...(x-  xn)  ^ 

(Xq  -  XjXXq  -  X^).  •  .  (  Xq  -  xj  0 
+  (x  -Xrt)(x  -X,).  . .  (x  -  x^) 

- U - i - a_  Y, 

(x  -Xf,)(x  -X  ).  . .  (x  -X  ,)  n 

n  u  n  1  n  n- 1 

which  yield  a  polynomial  of  degree  such  that  (x^)  =  Yq,  (x^)  =  y^^, 
. . . ,  V  (x^)  =  Y^,  as  is  easily  demonstrated  by  an  induction 

proof. 

When  X  is  within  the  limits  of  the  table,  the  points  x^,  x^ . x^ 

are  spaced  about  it  to  best  advantage,  but  if  x  is  outside  the  limits 
of  the  table,  the  nearest  n-t  1  points  in  the  table  are  used  in  the 
formula. 


I 
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g.  RDF. 


TITLE: 

ID: 

Classification: 

Programers: 


Read  Data,  Fixed  Formats 
RDF 

Input  Routine 

W.  Silverman  and  R,  £.  Mann 


PURPOSE; 

To  load  data  from  Hollerith  cards  through  the  088  card 
reader,  or  from  80  character  records  on  BCD  tape.  The  data 
are  on  symbolic  cards  of  the  types  used  for  the  pseudo  operations 
DEC,  OCT,  and  BCD  as  described  In  the  assembly  routine. 


USAGE: 


1.  Calling  Sequence 
ENA  A 
ENQ  B 

a  RTJ  RDF 

a+1  NOP  LiLTN) 

a+Z  Error  Return 

a  +  3  Normal  Return 

L(LTN)  =  Location  of  Logical  Tape  Number 
LTN  (Logical  Tape  Number)  =  1-48  for  tape  Input 

=  49  for  card  input 


After  a  description  of  the  input  formats,  the  function  of  A  and  B 
will  be  discussed. 


2.  Input  Formats: 

Load  Is  controlled  by  a  symbolic  operation  code  in  columns 
10, 11,  IZ  of  the  Input  cards  (or  records).  There  are  five  permis¬ 
sible  operation  codes. 


1.  SLJ--This  operation  always  causes  loading  to  be  termin¬ 
ated.  A  transfer  address  may  appear  In  decimal  or  octal 
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beginning  In  column  20.  The  address  is  assumed  decimal 
if  it  is  followed  by  a  D  or  blank,  and  octal  if  it  is  followed 
by  a  B.  If  index  modification  of  the  transfer  address  is 
desired,  column  17  or  18  may  contain  an  index  register 
number.  Control  is  transferred  to  the  (modified)  transfer 
address,  or  if  there  is  no  transfer  address,  control  is 
returned  to  a  t  3  in  the  calling  sequence. 

2.  REM--  This  record  will  be  skipped.  No  conversion  or 
operation  results.  REM  cards  may  be  used  as  spacers  or 
tags  in  the  data  deck. 

Conversion  operations: 

3.  BCD — n  words  of  binary  coded  decimal  information 
beginning  in  column  21  are  loaded  into  consecutive  memory 
locations,  n  is  in  column  20,  lCn^7.  The  information 
runs  through  column  20  -t-  8n. 

4.  OCT--Octal  data  beginning  in  column  20  and  terminating 
with  the  first  blank  column  are  interpreted  as  octal  integers, 
and  converted  to  binary  integers.  Successive  words  are 
separated  by  commas  and  loaded  into  consecutive  locations. 
Each  word  may  consist  of  a  -t-  or  -  sign  and  up  to  16  octal 
digits.  If  no  sign  appears,  a  ■¥  sign  is  assumed. 

5.  DEC --Decimal  data  beginning  in  column  20  and  termin¬ 
ating  with  the  first  blank  column  are  loaded.  Successive 
words  are  separated  by  commas  and  loaded  into  consecutive 
locations.  Each  word  may  consist  of  a  sign,  -t-  or  -  or  none, 
up  to  15  decimal  digits  with  a  decimal  point  if  desired,  a  D 
or  E  followed  by  a  signed  or  unsigned  decimal  scale  factor, 
and  a  B  followed  by  a  signed  or  unsigned  binary  scale 
factor.  Presence  of  a  binary  scale  factor  will  cause  the 
number  to  be  loaded  as  a  fixed  point  decimal  number  with  the 
binary  point  to  the  right  of  the  bit  position  given  by  the  binary 
scaling.  If  a  decimal  scale  factor  or  decimal  point  is  present, 
the  number  will  be  loaded  as  a  fixed  point  number  as  above  if 
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a  binary  scaling  is  present,  or  as  a  (scaled)  floating  point 
number  if  no  binary  scaling  is  present.  If  a  binary  scaling 
is  present,  it  must  lie  between  -47  and  47,  and  the  decimal 
scaling,  if  any,  must  lie  between  -28  and  28.  If  no  decimal 
point  or  scale  factor  is  present,  the  number  will  be  loaded 
as  an  integer.  A  few  illustrative  examples  follow: 

1)  1.2345  will  be  loaded  as  floating  point  1.2345 

2)  12345E-3  will  be  loaded  as  floating  point  12.  345 

3)  1.  2345D-2  will  be  loaded  as  floating  point  .  012345 

4)  12345  will  be  loaded  as  integer  12345 

5)  12345B15  will  be  loaded  as  fixed  point  12345.  0  with 
binary  point  to  the  right  of  bit  position  15, 

6)  12.  345D2B13  will  be  loaded  as  fixed  point  1234.  5 
with  binary  point  to  the  right  of  bit  position  13. 

7)  12345B-3  will  be  loaded  as  fixed  point  12345,0  with 
binary  point  to  the  right  of  bit  position  -3,  that  is  as 
(rounded)  integer  1543  =  12345/8  to  the  nearest  integer. 

Any  data  card  BCD,  OCT,  or  DEC  may  have  an  absolute 
numerical  address  in  columns  1-8.  The  data  on  the  card  is 
loaded  relative  to  that  address.  The  address  is  treated  as 
octal  or  decimal  according  to  the  same  conventions  used  for 
the  transfer  address  on  the  SLJ  card, 
e.  g.  The  card 

2000B  BCD  2ABCDEFGHIJKLMNOP 

will  be  loaded  so  that: 

(2000  B)  =  ABCDEFGH 
(2001B)  =  IJKLMNOP 

A  and  B  (the  addresses  in  the  accumulator  and  quotient 
registers  upon  entry  to  RDF)  condition  where  the  data  on 
cards  with  blank  address  fields  is  loaded  and  the  number  of 
words  actually  loaded  by  RDF. 

If  A  /  0,  the  first  card  with  a  blank  address  field  is  loaded 
relative  to  A,  the  data  going  into  addresses  A,  A-f  1, . . . , 
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A-t-n-1,  where  n  is  the  number  of  data  words  on  the  card.  The 
data  from  the  next  card  with  a  blank  address  field  is  loaded 
relative  to  A  into  locations  A-f  n,  A-t  n-f  1, . . . ,  A-fnfm-1 
where  m  is  the  number  of  data  words  on  the  card.  And  so  on, 
each  card  with  a  blank  address  field  being  loaded  relative  to  A. 

If  A  =  0,  the  first  card  must  have  an  absolute  numerical 

address  in  columns  1>8.  Loading  proceeds  relative  to  the  last 
such  address.  Thus  if  the  first  data  card  has  the  address  2000B, 
the  second  card  has  a  blank  address  field,  the  third  has  the 
address  3000B,  and  the  remaining  data  cards  have  blank  address 
fields,  the  data  from  the  first  cards  will  go  into  consecutive 
address'  s  starting  at  3000B.  If  B  ^  0,  B-A-f  1  words  will  be 
loaded  unless  the  load  is  terminated  by  an  SLJ  card  or  an  error 
(see  errors  below). 

If  B  s  0,  an  indefinite  number  of  words  will  be  loaded  until 
the  routine  is  terminated  by  an  SLJ  card  or  an  error. 

3.  RDF  requires  397  locations 

4.  RDF  uses  19  conrunon  erasable  storage  locations 

5.  Errors: 

1.  End  of  file  on  tape  or  card. 

The  end  of  file  flag  (77713B)  is  set  non-zero  and  if  no 
other  errors  occurred,  control  is  returned  to  a  -f  3 
in  the  calling  sequence. 


2.  Parity  error  on  tape. 

The  RTT  flag  (77712B)  is  set  non-zero,  and  the  routine 
continues,  although  the  conversion  of  the  record  on  which 
the  parity  error  occurred  may  be  inaccurate.  When  the 
routine  is  terminated,  control  is  returned  to  the  error 
return,  a  4  2,  in  the  calling  sequence,  with  the 
number  of  records  which  had  parity  errors  in  the  upper 
address  of  the  A  register,  and  the  transfer  address,  if 
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any,  In  the  Q  register. 

3.  Illegal  punch  on  card. 

The  RTT  flag  is  set  with  I's  in  bits  0  through  39  corres¬ 
ponding  to  erroneous  columns  on  the  card  (these  may  apply 
to  the  first  or  second  half  of  the  card,  but  not  both).  Control 
is  returned  to  a  -f  2  in  the  calling  sequence,  with  the  upper 
address  of  A  set  to  1. 

4.  Format  error: 

This  may  be  caused  by  an  illegal  operation  code  in 
columns  10,  11,  12;  an  illegal  address;  an  illegal  character 
in  column  17  or  18  of  an  SLJ  card;  or  an  illegal  character  in 
column  20  of  a  BCD  card.  Control  is  restored  to  a  -t-  2  in 
the  calling  sequence  with  bit  47  of  A  set  to  1 . 

5.  Data  error: 

This  may  be  caused  by  an  illegal  character  in  a  numerical 
field;  more  than  16  digits  in  an  octal  numerical  field,  or  more 
than  14  in  a  decimal  numerical  field;  more  than  one  sign  in  a 
numerical  field  or  a  sign  which  is  preceded  by  a  number  in 
the  field;  more  than  one  decimal  point  in  a  decimal  numerical 
field;  too  large  a  scale  factor;  a  decimal  number  which  when 
scaled  does  not  fit  the  A  register,  or  which  is  too  large  or 
too  small  for  the  floating  point  format 

,1023  .  ,  1023k 

(^-2  or  ^  -2-  ). 

In  case  of  a  data  error,  the  word  in  storage  correspond¬ 
ing  to  the  erroneous  field  is  set  equal  to  minus  zero,  and 
loading  continues.  At  the  end  of  the  routine,  control  is  re¬ 
turned  to  a  -t-  2  in  the  calling  sequence,  with  the  number  of 
erroneous  fields  in  the  lower  address  of  the  A  register. 

6.  Termination  of  loading: 

An  HLJ  card,  an  end-of-file,  or  a  format  error  auto¬ 
matically  terminates  loading.  If  B  0  (in  the  Q  register 
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of  entry),  loading  is  terminated  when  B-A-t-  1  words  have 
been  loaded.  The  program  then  hunts  forward  through  the 
remaining  cards  (or  records  on  tape)  to  the  first  SLJ  card, 
or  end-of-file.  The  tape  or  card  reader  is  left  set  to  read 
the  next  succeeding  record  or  card.  Normal  return  in  this 
case  is  to  a  -(■  3  in  the  calling  sequence. 

7.  All  conversions  are  performed  rapidly  enough  so  that 
the  1607  tape  mechanism  (or  the  088  card  reader)  can  be 
operated  at  full  speed;  3000  records  per  minute  on  tape  or 
650  cards  per  minute  tl .rough  the  088. 

8.  Floating  point  numbers  with  large  negative  exponents  may 
be  accurate  to  only  35  bits. 
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h.  EfltAt£ 


MATRIX  ROTATION 


The  following  transformation  was  used  as  a  subroutine  to  transform 
from  one  Cartesian  coordinate  system  to  another  Cartesian  coordinate  system. 


Z,  Z' 


I .  Rotation  of  the  initial 
system  by  an  angle  p 
counterclockwise  about 
the  Z  axis. 


2,  Rotation  of  the  primed 
system  by  the  angle  6 
counterclockwise  about 
the  X'  axis. 


3.  Rotation  of  the  double  primed 
system  by  the  angle  4  counter¬ 
clockwise  about  the  Z"  axis. 


A 


cos  4  cos  f  -  cos  6  sin  9  sin  4 
cos  4  vln  9  +  cos  6  cos  9  sin  4 
sin  6  sin  4 


-sin  4  cos  9 -cos  6  sin  9  cos  4  sin  6  sin  9 
-  sin  4  sin  9  -fcos  6  cos  9  cos  4  -  *in  6 cos  9 
sin  6  cos  4  cos  6 
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This  program  is  callable  by  FORTRAN. 
To  use: 

CALL  ROTATE  (A,  B.  C,  D,  E) 


where 


A 

B 

C 

D 

E 


The  angle  7 
The  angle  6 
The  angle  ^ 

The  input  3  component  vector 
The  output  3  component  vector 


I 
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(i 

i.  s£xi:a& 

The  Subroutine  SETTAB  is  used  to  set  up  a  table  of  print  times. 
This  table  includes  both  even  increment  print  times  and  abnormal  print 
times  such  as  ignition,  drop  stages,  and  burnout. 

The  subroutine  is  callable  by  FORTRAN.  To  use: 

CALL  SETTAB  (NS,  TMI,  TMCC,  TMBO,  DT,  TABLE,  IX) 

where 


NS 

= 

number  of  stages 

TMI 

= 

table  of  ignition  times 

TMCC 

= 

table  of  coefficient  change  times 

TMBO 

= 

table  of  burnout  time 

DT 

= 

even  time  increments 

TABLE 

s 

resulting  table 

IX 

= 

size  of  the  table 
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j,  TWO  BOD 


SYMBOLS 


BNqLtSa 


a 


b 

C 


*E 

E 

£ 

CM  - 


- 


h 

z 


h 


i 

k 

M 


P 

R 

R 

R' 

S 


t 

o 


Semi- major  axis  of  orbit 

Equatorial  radius  of  the  earth 

Semi- minor  axis  of  orbit 

Earth  parameter  defined  by  equation  3.1 

Eccentricity  of  the  orbit 

2  2 

Eccentricity  of  the  earth  =  2f  -  f  ) 

Eccentric  anomaly 
Flattening  of  the  earth 
Gravitational  constant  of  the  earth 
Geodetic  altitude 

Angular  momentum  about  subscripted  axis 
Total  a-igular  momentum 
Apogee  altitude 
Petigee  altitude 

Energy  parameter  defined  by  equation  2-2 
Inclination  of  orbit  plane  to  equatorial  plane 
Constant  used  to  obtain  period  {Z%/ -JCM) 

Mean  anomaly 
Period  of  the  orbit 

Geocentric  radius  from  center  of  earth  to  vehicle 
Change  in  R  with  respect  to  time 
Geocentric  radius  used  for  earth's  velocity 
Earth  parameter  defined  by  equation  3-2 
Time  of  entry  into  orbit 


(NM) 

(NM) 

(NM) 

(NM) 


(RAD) 

(FT^/SEC^ 

(NM) 


(NM) 

(NM) 

(DEG) 


(MIN) 

(NM) 

(NM/SEC) 

(NM) 

(NM) 

(SEC) 
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*b‘  ^b- 

X,  y,  z 


x',  y',  z' 


X",  y",  z" 
X.  Y,  Z 


■GR£EIC 


P 

Y 


u 


e 

c 

I 

e 

c 


V- 

V 


SYMBOLS 

-  Time  of  perigee  passage 

-  Time  at  point  in  orbit 

-  Velocity  of  vehicle  with  respect  to  the  earth 

~  Velocity  of  earth's  surface  below  vehicle 

-  Velocity  of  vehicle  with  respect  to  non¬ 
rotating  earth 

-  Earth  centered  coordinates  of  the  space  vehicle 

-  Geodetic  coordinate  system  at  vehicle 

-  Geodetic  coordinate  system  below  vehicle  on 
surface  of  the  earth 

-  Coordinate  system  in  orbit  plane 

-  Earth  centered  coordinate  system 

-  Coordinates  of  the  look  station 


-  Angle  between  geodetic  east  axis  and  X  axis  of 
coordinate  system  tangent  to  a  geodetic  earth; 
also  the  aspect  angle 

-  Initial  azimuth  (positive  c.  w.  from  north) 

-  Initial  flight  path  angle  (positive  up) 

-  Inertial  flight  path  angle  (positive  up) 

-  Argument  of  perigee 

-  Rotational  rate  of  the  earth 

-  Longitude  of  ascending  node 

-  Longitude 

-  Geodetic  latitude 

-  Geocentric  latitude  from  vehicle 

-  Geocentric  latitude  from  earth's  surface 

-  Argument  of  latitude 

-  True  anomaly 


(SEC) 

(SEC) 

(EPS) 

(EPS) 

(EPS) 
(Eigure  3) 
(Eigure  4) 

(Eigure  4) 
(Eigure  Z) 
(Eigure  1) 
(Eigure  3) 


(DEG) 

(DEG) 

(DEG) 

(DEG) 

(DEG) 
(RAD/ SEC) 

(DEG) 

(DEG) 

(DEG) 

(DEG) 

(DEG) 

(DEG) 

(DEG) 
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SUBSCRIPTS 


b 

c 

E 

G 

i 

I 

o 

R 

s 

X.  Yt  as 
1 

Z 


- 

refers 

to 

- 

refers 

to 

- 

refers 

to 

- 

refers 

to 

- 

refers 

to 

- 

refers 

to 

- 

refers 

to 

refers 

to 

- 

refers 

to 

' 

refers 

to 

refers 

to 

plane 

- 

refers 

to 

space  vehicle 
geocentric 
the  earth 
geodetic 

point  in  orbit 
nonrotating  earth 
Initial  or  burnout  values 
rotational  velocity  radius 
look  station 

X,  Y,  Z,  coordinate  system 

coordinate  system  at  station  parallel  to  equatorial 
coordinate  system  at  look  station 


A  dot  over  a  variable  indicates  the  time  derivative  of  that  variable. 

A  delta  (A)  in  front  of  a  variable  indicates  a  difference  in  that  variable. 
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1.  INTRODUCTION. 

a.  Input ■. 

The  inputs  are  a  position  and  velocity  vector  in  a  right-hand  Earth- 
centered  coordinate  system  with  the  X  axis  along  the  node  of  the  equatorial 
plane  and  the  Greenwich  meridian  and  the  Z  axis  along  the  north  polar  axis. 
The  number  of  time  increments,  changes,  and  look-angle  stations  must  be 
in  COMMON.  To  use  the  subroutine,  control  number  4  must  be  set  equal 
to  1  or  2. 

b.  Printout  time  changes. 

Since  some  portions  of  the  trajectory  are  more  important  than 
others,  provisions  are  included  for  changing  the  time  increment,  with  the 
use  of  a  maximum  of  five  different  time  increments.  If  the  ellipse  intersects 
the  earth,  the  trajectory  computation  stops  at  this  time.  If  it  does  not  inter¬ 
sect  the  earth,  the  computation  will  continue  for  a  pre specified  number  of 
orbits. 

The  program  computes  the  classical  orbital  elements  (a,  b,  e,  P, 

*  ^p*^**^'  ''o’  ^o  *• 

After  the  orbital  elements  are  computed,  a  matrix- rotation  is  set  up 
to  transfer  from  the  orbital  plane  to  the  equatorial  plane  coordinate  system. 
Time  is  incremented  and  a  new  true  anomaly  is  obtained. 

The  corresponding  orbit  plane  coordinates  are  then  transformed  by 
a  matrix  rotation  to  the  equatorial  coordinate  system.  The  latitude  and 
longitude  are  obtained  by  taking  "arc  tangents"  and  adding  the  effects  of  a 
rotating  earth.  These  values  are  stored  along  with  the  radius  vector  to  the 
vehicle  and  are  used  to  find  the  "look  angles"  at  various  stations. 

The  matrix  subroutine  described  in  rotate  is  used  for  the  coordinate 
transformation. 
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(3) 

(4a) 

(4b) 

(4c) 


i  =  cos"^  (h  /h) 

z 

h 

-1  ( — S-) 

Q  =  tan  -h  It  (-b  )  is  negative, 

y  y 

Q  =  D  +  180° 


(4d) 

(5) 

(6) 

(7) 

(8) 
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cos 

V  = 

o 

RGM 

(9) 

sin 

V  = 

o 

hR 

GM 

(10) 

V 

o 

= 

_ ,  /  sin  V  \ 

tan"  1  - 2_  1  If  (cos  V  )is  negative, 

\COS  V  /  ° 

°  V  =  V  +  180° 

o  o 

(11) 

cos 

(Y  h,  -  X  h^, 

(12) 

h 

= 

tan"^  (  — — —  j  If  (cos  n)  is  negative 
Vcos  ji  / 

'  '  (i  =  +  180° 

(13) 

u 

= 

(14) 

b 

= 

/  2 
a  vl  -  e 

(15) 

P 

= 

(16) 

‘^A 

= 

a  (1  +  e)  -  aj. 

(17) 

= 

a  (1  -  e)  -  a^. 

(18) 

E 

o 

= 

2  tan"  ^  '^) 

(19) 

t 

T  = 

(£  -  e  slnE  )P 
o  o 

(20) 

o 

2k 
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d.  Trajectory  points. 

M  =  -^  (tj  -  T)  <22) 


E 


2 


(E  j  -  e  sin  E  j  -  M) 
1  -  e  cos  Ej 


Iterate  until 


V 


i 


2  tan' 


-1 


tan 


5 

2 


) 


(23) 


(24) 


B  -  al  l--:.  S  l _ 

i  1  +  e  cos  V. 


Vi 


I  pe  sin  Vj 

n  =  TTTToTT^ 


X  =  cos  V. 


y  =  sin  V. 


a"  =  0 


(25) 

(26) 

(27) 

(28a) 

(28b) 

(28c) 
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X 

Y 

Z 


cos  0  -  sin  Q 
sin  Q  cos  Q 
0  0 


0 

0 

1 


1 

0 

0 


0  0 
cos  i  -sin  i 
sin  i  cos  i 


cosu  -sinu 
sinui  cosu 
0  0 


II 

X 

y 

II 

Z 

(29) 


X  =  x"  (cos  Q  cosu  -  sin  Q  cos  i  sinu)  -  y"  (cos  Q  sinu  sin  Q  cos  i  cosu) 
Y  =  x" (sin  Q  cosu  +  cos  Q  cos i  sinu)  -  y"  (sin  Q  sinu  -  cos  Q  cos  i  cosu) 
Z  =  x"  (sin  i  sinu)  +  y"  (sin  i  cosu) 

(30) 


6 

c 


tan 


(31) 


If  X  is  negative,  9^.  =  ^  180“ 


(32) 


’c  =  ’c  -  “E  <‘i  -  ‘o’ 


(33) 


e.  Geocentric  to  geodetic. 


f.  Range. 


The  great  circle  range  is  obtained  by  computing  the  range  angle.  This 
angle  is  obtained  by  taking  the  dot  product  of  the  launch  vector  with  the  radius 
vector  at  any  given  time. 


R.A.  =  Cos'^  {(X^  X  +  Yj^  Y  +  Z)/  v^X^  +  yJ^  +  Z^  v^X^  +  Y^  +  Z^} 

Range  *  Rj,  (avg)  (R.A. ) 
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8*  ItOOK  >ngUgi 

SiaiiSOL  given  station  (h^  ,  6^, ,  «p^ ) 


S  8 


C  = 


*  (1  -  sin^e^) 


2  .4  2c  *Vi  ‘^s-^'s'^-'E 


S_  =  C(1  -  e  ‘  ) 


(37) 


|q-[.vv 


)  sin  6^  -t-  (C  -I-  h  )  cos  6- 

G  s  S  G 

S  *8 


.]" 


(38) 


(39) 


®c 


<y.  =  0 


z  =  sin  6 
H  E  c 


(40) 


Yk  = 


R^  cos  6^  cos  A  cp 
R.  cos  6  sin  A  q> 

b  Cfa  ^ 

R,  sin  0 


where  Atp  =  «)>(,-«?, 


(41) 


R.  =  Rr  +  P 


P  =  \  -  \ 


(42) 
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P 


yi 

*1 


cof  9  cos  A®  -  R_  cos  6_ 
%  ^  E  c 


cos  6  sin  A® 

% 


sin  8  -  R„  sin  6 

c.  c 


(43) 


*2 

cos  6q  0  -  sin  0^ 

*1 

yz 

0  1  0 

yi 

*2 

sin  0Q  0  cos  6q 

*1 

Elevation  ang.  =  tan 


Azimuth  ang.  =  tan 


Ifx^neg. 
add  180° 


If  (azimuth)  is  neg,  add  360° 


(44) 


(45) 


(46) 


h.  Aspect  angle. 

The  aspect  angle  is  defined  as  the  angle  between  the  spin  axis  of  the 
space  vehicle  and  the  vector  from  the  look  station  to  the  vehicle.  This  program 
assumes  the  vehicle  remains  fixed  in  space.  The  angle  is  obtained  by  the 
following  equations: 

a  =  Cos'^jlXXj  +  Yy^  +  Zz^)/  y/i}  +  /  y/x^  +  yj  +  z^ 

(47) 


1Z9 


VARIABLES  USED  IN  TWO  BODY  SUBROUTINE 
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ALT 

Geocentric  radius  to  vehicle 

AMU 

Argument  of  latitude 

ANS 

Output  vector 

ANT 

Aspect  angle 

APOGEE 

Apogee  altitude  of  the  orbit 

ARA 

Average  radius  of  the  Earth 

AX 

X  Component  of  inertial  burnout  position  vector 

AXDOT 

X  Component  of  inertial  burnout  velocity  vector 

AY 

Y  Component  of  inertial  burnout  position  vector 

AYDOT 

Y  Component  of  inertial  burnout  velocity  vector 

AZ 

Z  Component  of  inertial  burnout  position  vector 

AZDOT 

Z  Component  of  inertial  burnout  velocity  vector 

AZIM 

Azimuth  of  look  vector 

A1 

Burnout  position  vector 

A2 

Burnout  velocity  vector 

AlA 

Not  used  in  two  body  ♦ 

AIB 

Not  used  in  two  body  * 

C 

^  Ree/  %/ \  -  e^  sin  6^ 

CAPO 

Longitude  of  ascending  node 

CLAL 

Cosine  of  the  launch  latitude  * 

CLOL 

Cosine  of  the  launch  longitude 

CMU 

Cosine  of  MU 

CSLANT 

Cosine  of  inclination  angle 

CT 

Cosine  of  geocentric  latitude 

CTRVA 

Cosine  of  true  anomaly 

DEL TAT 

Stored  time  increment  * 

DPHI 

Difference  in  longitude 

DUMTB 

Not  used  in  two  body  * 

E 

Eccentricity  of  the  Earth 

*  Stored  in  COMMON 
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ELEV 

Elevation  of  look  angle  vector 

ENDT 

End  of  time  increment  * 

ER 

Size  of  allowable  error 

ETA 

Anomaly  =  period/ 2n 

EX 

Orbital  eccentricity 

EXANOM 

Eccentric  anomaly 

EXIMP 

Impact  eccentric  anomaly 

El 

Used  to  compute  eccentric  anomaly 

E2 

Used  to  compute  eccentric  anomaly 

F 

Flattening  of  the  Earth 

FMIN 

Output  -  time  in  minutes 

FMINF 

Function  for  obtaining  FMIN 

FPA 

Flight  path  angle 

FX 

Output  vector  from  Rotate 

CM 

Earth  gravitational  attraction  constant 

GMl 

Earth  gravitational  attraction  constant 

H 

Angular  momentum  parameter 

HH 

Angular  momentum  squared 

HOUR 

Output  time  in  hours 

HOURH 

Function  for  obtaining  HOUR 

HX 

X  Component  of  angular  momentum 

HY 

Y  Component  of  angular  momentum 

HZ 

Z  Component  of  angular  momentum 

HI 

Angular  momentum 

I 

Utility  index 

IT 

Utility  index 

J 

Utility  index 

K 

Utility  index 

KBATT 

Number  of  Batt  output  tape  * 

MEX 

Utility  index 

MM 

Utility  index 

NAME 

Name  stored  in  common  * 

*  Stored  in  COMMON 
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NAT 

NOST 

NOT 

NIC 

ORBT 

PERIGEE 

PERIOD 

pm 

pmx 

PI 
PN 
POP 
PR  ALT 
PRLAT 
PRLON 
PR  TIME 
R 

RAD 

RANGE 

RDOT 

RE 

REE 

RF 

RLAU 

RS 

S 

SEC 

SECF 

SHT 

SLAL 

SLANT 

SLOL 

SLR 


Number  of  output  tape  * 

Number  of  look  angle  etatlone  * 

Number  of  time  changes  * 

Not  used  In  two  body  * 

Number  of  orbits  of  orbital  vehicle  * 

Perigee  altitude  of  the  orbit 

Time  to  make  one  revolution  of  orbit 

Stored  longitude 

Used  to  compute  longitude 

=  It  =  3.141592654 

Page  count 

Earth  flattening  constant 
Geodetic  latitude  of  the  vehicle 
Geodetic  altitude  of  the  vehicle 
Longitude  of  the  vehicle 
Time  from  launch 
Length  of  position  vector  n.  m. 

=  n  /ISO.  =  1/57.2957795 
Great  circle  range  on  surface  of  Earth 
Rate  of  change  along  R  with  respect  to  time 
Radius  of  Earth  =  f(6c) 

Equatorial  radius  of  Earth  -  n.  m. 

Length  of  position  vector  (ft) 

Radius  to  launch  point 
Radius  of  look  station 
AC  (1  -  e^) 

Output  time  In  seconds 

Function  for  obtaining  output  time  in  seconds 

Look  station  altitude 

Sine  of  the  launch  latitude  * 

Orbital  inclination  angle 

Sine  of  the  launch  longitude 

Distance  from  tracking  station  to  vehicle 


*  Stored  in  COMMON 
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SMAJ 

Semi  major  axle 

SMIN 

Semi  minor  axis 

SMON 

Argument  of  perigee 

SPH 

Look  station  longitude 

SSLANT 

Sine  of  inclination  angle 

ST 

Sine  of  theta 

STH 

Look  atation  latitude 

STHC 

Look  station  geocentric  latitude 

STRVA 

Sine  of  the  true  anomaly 

SX 

Orbital  plane  coordinate  system  vector 

SI 

Variables  used  to  compute  output  variables 

S2 

Variables  used  to  compute  output  variables 

S3 

Variables  used  to  compute  output  variables 

S4 

Variables  used  to  compute  output  variables 

S5 

Variables  used  to  compute  output  variables 

SP6 

Used  to  calculate  aspect  angle 

TAG 

Not  used  in  two  body  * 

THETA 

Geocentric  latitude 

TIME 

Independent  variable 

TIMP 

Time  to  impact 

TINCR 

Time  increment 

TOPP 

Time  of  perigee  passage 

TRANOM 

Initial  true  anomaly 

TRIMP 

Impact  point  true  anomaly 

trunu 

True  anomaly 

VELOC 

Velocity  at  any  point  of  the  orbit 

VI 

Initial  inertial  velocity 

WE 

Rotational  velocity  of  the  Earth 

X 

X  Component  of  look  vector 

XAZ 

Azimuth  of  look  vector 

XEL 

Elevation  of  look  vector 

XIM 

X  Component  of  position  vector  for  range  computations 

*  Stored  In  COMMON 
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XLAV 

X  Component  of  launch  vector 

XLLO 

Launch  longitude  * 

XM 

Mean  anomaly 

XP 

Used  in  matrix  rotation 

XSHT 

Stored  tracking  station  altitude  * 

XSPH 

Stored  tracking  station  longitude 

XSTH 

Stored  tracking  station  latitude  * 

XX 

=  y/ ^  ,  storage  variable 

XXI 

Sine  of  station  geodetic  latitude 

XX2 

Cos  of  Station  geodetic  Latitude 

XX3 

Used  in  station  radius  calculation 

XX4 

Used  in  station  radius  calculation 

Y 

Y  Component  of  look  vector 

YIM 

Y  Component  of  position  vector  f( 

YLAU 

Y  Component  of  launch  vector 

YY 

Storage  variable 

Z 

Z  Component  of  look  vector 

ZIM 

Z  Component  of  position  vector  f( 

ZIZ 

Range  angle 

ZLAU 

Z  Component  of  launch  vector 

ZP 

Used  in  matrix  rotation 

Z2 

Number  of  lines  per  page  * 

*  Stored  in  COMMON 
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k.  WOTF.  SAUF. 

A.  IDENTIFICATION 


TITLE; 
CO-OP  ID: 
CATEGORY: 
Programer: 
Date: 


WRITE  ON  TAPE 
WOTF 

General  Tape  Handler 
J.  W.  Vise 
August  4,  1961 


B.  PURPOSE; 

This  routine  Is  used  to  write  a  record  of  arbitrary  length  on  a 
magnetic  tape  in  either  binary  or  BCD  mode 


C.  USAGE: 

1.  Calling  Sequence 

CALL  WOTF  (a.  b,  n,  m) 


n 

s 

1-48,  MT 

m 

0,  BINARY 

m 

= 

1,  BCD 

a 

= 

First  LOCN 

b 

= 

Last  LOCN 
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CLASSIFICATION:  Utility  Routine 
TITLE:  SAVE  TAPE 
ID:  SAVF 

PURPOSE:  To  Inform  the  operator  that  a  tape  is  to  be  saved  and  to 
rewind  the  tape  with  Interlock.  E.  O.  F.  and  98ssssssENDs8sss 
0Fs88888TAPE8888  l8  Written  on  the  tape.  (Here  8  repre8ent8 
blanks). 

USAGE: 

1.  Calling  Sequence: 

This  routine  Is  used  with  the  following  FORTRAN  calling 
sequences: 

CALL  SAVF  (N),  or 
DUMMY  =  SAVF  (N) 

where  N  is  the  number,  1-48,  of  the  logical  tape  to  be  saved. 
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7.  PROGRAM  LISTING 
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***  SPURT  *** 

SEPTEMBER  19,  1962  MOO  OF  AUG.  21,  1962 
IF  KNTRL  (l)“l,  USE  METRO  DATA 
IF  KNTRL  (l)-2,  USE  WINDS  ONLY 
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SUBROUTINE  IMPACT(jj) 

THIS  SUBROUTINE  COMPUTES  THE  TRAJECTORY  OF  THE  EMPTY  STAGES 
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SP.UBT  SAMPLE.P.RlMiaiLX.PAIA 


1.  INPUT  DATA  -  Launch  Parameters  and  Spin  Table 
Z.  INPUT  DATA  -  Stage  Weight  and  Aerodynamic  Parameters 

1  page  per  Stage 

3.  SPURT  OUTPUT  DATA 

4.  SPURT  OUTPUT  DATA 

5.  SPURT  OUTPUT  DATA 

6.  SPURT  OUTPUT  DATA 

7.  TWO-BODY  OUTPUT  DATA  -  Orbital  Elements 

8.  TWO-BODY  OUTPUT  DATA  -  Keplerian  Trajectory 

9.  TWO-BODY  OUTPUT  DATA  -  Look-Angles 

10.  IMPACT  OUTPUT  DATA  -  Position  Information 

11.  IMPACT  OUTPUT  DATA  -  Velocity  Information 
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DEFINITIONS 


DEFINITIONS;  All  trajectories  ore  with  respect  to 

an  oblate  earth. 

RANGE  COORDINATE  SYSTEM;  A  left  hand  orthogonal  cartesian 

coordinate  system  with  the  XY 
plane  tangent  to  an  oblate  earth  at 
the  launch  point.  The  Z  axis  is 
positive  upward  along  the  local 
vertical  and  the  X  axis  is  positive 
along  the  launch  azimuth. 


Laundt 

Asiauth 


TIME  (Sec): 


LATITUDE  (Deg): 


LONGITUDE  (Dsf): 


Relates  to  the  first  motion  of  the 
vehicle  off  the  launch  site. 

Angle  between  the  normal  to  the 
reference  spheroid  passing  through 
the  vehicle  and  the  equatorial  plane. 
Positive  in  the  northern  hemisphere. 

Angular  distance  measured  from  the 
loot  of  the  Greenwich  meridian  to  the 
vehicle  sub-point  meridian.  West  of 
Greenwich  is  negative  and  East  is 
positive. 
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ALTITUDE 

VELOCITY  dpi); 

AZIMUTH  tC«g}; 

F..F..A.  tPtg): 

RANGE 

PEfLECTIOM.tW.MJ: 
PM  (Dgf); 

THETA  (Dey); 

THRUST  (Lhi): 

WEIGHT  (Lba): 


Distance  from  vehicle  to  the  surface 
of  a  geodetic  earth  below  the  vehicle. 

Velocity  of  vehicle  with  respect  to 
a  point  on  a  rotating  earth  directly 
below  the  vehicle. 

Angle  from  North  the  velocity  vector 
makes  with  the  local  meridian. 
Measured  positive  C.  W.  from  North. 

Flight  path  angle.  The  angle  the 
velocity  vector  makes  with  the  local 
horizon;  positive  upward. 

Arc  distance  from  launch  to  the  point 
under  the  vehicle  measured  along  the 
surface  of  the  earth. 

Distance  the  vehicle  has  deviated 
from  the  launch  axis  in  the  XY  plane. 

Euler  angle  between  the  longitudinal 
axis  of  the  vehicle  and  the  Z  axis 
in  the  range  coordinate  system 

Euler  angle  between  the  longitudinal 
axis  of  the  vehicle  and  the  XZ  plane 
measured  in  the  rotated  X'Y*  plane. 

Thrust  of  the  rocket  motor  with 
thrust  increase  due  to  pressure  de¬ 
crease  included. 

Weight  of  the  remaining  portion  of 
the  veMcle.  Given  in  terms  of  sea 
level  pounds. 
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DRAG 
MACH  NO; 


PERIOD  Ip): 

APOGEE  ALTITUDE  (h.): 


Absolute  total  acceleration  of  the 


vehicle  normalized  by  the  gravita¬ 
tional  constant  . 

o 


The  classical  dynamic  pressure  on 
the  vehicle  given  by  Si  of  the  density 
times  the  velocity  squared. 


Aerodynamic  drag  on  the  vehicle. 


Mach  number  of  the  vehicle. 


Position  vector  of  vehicle  in  range 
coordinate  system 

Velocity  vector  of  vehicle  in  range 
coordinate  system 


(Ref  3)  The  parameters  of  a  classi¬ 
cal  Keplerian  ellipse  based  on  an 


inverse  square  gravitational  field. 


One  half  of  the  longest  diameter  of 
the  Keplerian  ellipse. 


One  half  of  the  smallest  diameter  of  the 
Keplerian  ellipse. 


A  measure  of  the  flattening  of  the 
Keplerian  ellipse. 


The  time  that  a  space  vehicle  takes 
to  make  one  complete  orbit. 


The  distance  to  the  highest  point  on 
the  ellipse  from  the  surface  of  the 
Earth. 


The  radius  of  the  point  on  the  ellipse 
closest  to  the  Earth,  minus  the  radius 
of  the  Earth. 
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The  angle  between  the  orbit  plane 
and  the  equatorial  plane. 

The  angular  distance  measured  in 
the  orbit  plane  from  the  line  of 
nodes  to  the  line  of  apsides. 

The  angular  distance  measured  at 
burnout  time  from  Greenwich  east¬ 
ward  in  the  equatorial  plane  to  the 
point  of  intersection  of  the  orbit 
plane  where  the  vehicle  crosses 
from  south  to  north. 

The  angle  measured  at  burnout  time 
at  the  center  of  the  Earth  between 
the  line  of  apsides  and  the  radius 
vector  to  the  vehicle  measured  from 
perigee  in  the  direction  of  motion. 

The  angle  at  the  center  of  the  ellipse 
between  the  line  of  apsides  and  radius 
vector  of  the  auxiliary  circle  through 
a  point  which  has  a  projection  of  the 
ellipse  corresponding  to  the  initial 
true  anomaly. 

Time  in  hours,  minutes,  and  seconds 
of  the  vehicle  from  launch. 

The  velocity  of  the  vehicle  with 

f 

respect  to  a  coordinate  system  fixed 
to,  but  not  rotating  with,  the  Earth. 
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INERTIAL  FLIGHT  PATH  ANGLIL 

II2£S): 

LOOK  ANGLES; 

AZMUTH  (Peg); 

ELEVATION 

SLANT  RANGE  (N.M.); 

ASPECT  ANGLE  (Deg); 


The  angle  at  any  given  time  the 
inertial  velocity  vector  makes  with 
respect  to  the  perpendicular  to  the 
radius  vector  at  that  time. 

The  direction  to  position  a  tracking 
antenna  at  a  given  station. 

The  angle  that  the  projection  in  the 
horisontal  plane  of  the  vector  point¬ 
ing  to  the  vehicle  makes  with  the 
North  direction. 

The  angle  from  the  horiaon  to  the 
vector  pointing  to  the  vehicle. 

The  distance  from  the  tracking 
station  to  the  vehicle. 

Angle  between  the  vehicle  spin  axis 
and  a  vector  from  a  given  station  to 
the  vehicle. 
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